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ABSTRACT 

It is now generally agreed that multidimensional, multigroup, radiation hydrodynamics is an indis- 
pensable element of any realistic model of stellar-core collapse, core-collapse supernovae, and proto- 
neutron star instabilities. We have developed a new, two-dimensional, multigroup algorithm that can 
model neutrino-radiation-hydrodynamic flows in core-collapse supernovae. Our algorithm uses an ap- 
proach that is similar to the ZEUS family of algorithms, originally developed by Stone and Norman. 
However, we extend that previous work in three significant ways: First, we incorporate multispecies, 
multigroup, radiation hydrodynamics in a flux-limited-diffusion approximation. Our approach is capa- 
ble of modeling pair-coupled neutrino-radiation hydrodynamics, and includes effects of Pauli blocking 
in the collision integrals. Blocking gives rise to nonlinearities in the discretized radiation-transport equa- 
tions, which we evolve implicitly in time. We employ parallelized Newton-Krylov methods to obtain 
a solution of these nonlinear, implicit equations. Our second major extension to the ZEUS algorithm 
is inclusion of an electron conservation equation, which describes evolution of electron-number density 
in the hydrodynamic flow. This permits following the effects of deleptonization in a stellar core. In 
our third extension, we have modified the hydrodynamics algorithm to accommodate realistic, complex 
equations of state, including those having non-convex behavior. In this paper, we present a description of 
our complete algorithm, giving detail sufficient to allow others to implement, reproduce, and extend our 
work. Finite-differencing details are presented in appendices. We also discuss implementation of this 
algorithm on state-of-the-art, parallel-computing architectures. Finally, we present results of verification 
tests that demonstrate the numerical accuracy of this algorithm on diverse hydrodynamic, gravitational, 
radiation-transport, and radiation-hydrodynamic sample problems. We believe our methods to be of gen- 
eral use in a variety of model settings where radiation transport or radiation hydrodynamics is important. 
Extension of this work to three spatial dimensions is straightforward. 

Subject headings: hydrodynamics — radiative transfer — methods: numerical — supernovae: 
general — stars: collapsed — stars: interior 
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1. Introduction 

Development of a numerical description of neutrino-radiation-hydrodynamic phenomena presents nu- 
merous modeling challenges. These challenges result from the complexity of the system — complexity stem- 
ming both from the sheer number of important physical components, and from the diversity and variability 
of microphysical interactions among these components. In the core-collapse supernova problem, this com- 
plexity presents itself in a number of ways. There is complexity associated with the interacting flows of 
matter and neutrino radiation. These flows exhibit coupling strengths that vary spatially, and change as a 
region evolves. There is typically strong coupling between matter and neutrino-radiation in dense regions 
near the center, weak coupling in the more diffuse outer regions of the core, and some kind of intermediate- 
strength coupling elsewhere. Where coupling is weaker, neutrino radiation is not in local thermodynamic 
equilibrium (LTE) with the matter. Adding to this is the presence of multiple species of neutrinos that are 
coupled through pair-production processes and through the exchange of energy and lepton number with 
matter. Hence, neutrino distributions span a large range of classical and quantum-mechanical behavior. 
There is also a complex nuclear chemistry, mediated by both the strong and weak interactions, present in the 
problem. Finally, the matter is described by a non-ideal-gas equation of state (EOS), which includes phase 
transitions and other complex behavior. 

Modeling such a phenomenon requires a set of neutrino-radiation-hydrodynamic equations that ac- 
counts for all the aforementioned complexities, which must then be integrated forward in time from an 
initial model. In this paper, we present an algorithm for accomplishing this in two spatial dimensions (2-D). 
The extension of this work to three spatial dimensions (3-D) is straightforward, and long-timescale 3-D sim- 
ulations of core-collapse supernovae using this approach should be computationally tractable in the next few 
years. In the sections that follow, we provide a description of our algorithm, supplying detail sufficient for 
others to implement and replicate this work. Although other algorithms for solving 1-D neutrino-radiation- 
hydrodynamics equations, or portions of these equations, have been published (see, for example, Yueh & 
Buchler 1977; Schinder & Shapiro 1982; Bruenn 1985; Myra et al. 1987; Schinder 1988; Schinder, Blud- 
man, & Piran 1988; Schinder & Bludman 1989; Mezzacappa & Bruenn 1993b; Swesty 1995), the only, al- 
beit partial, published descriptions of an algorithm for solving multidimensional Eulerian neutrino-radiation 
equations are those of Miller, Wilson, & Mayle (1993) and Buras et al. (2006). 

In our development of a radiation-hydrodynamic algorithm relevant to supernovae, we have drawn 
heavily from work originally performed for the development of the ZEUS family of multidimensional 
Eulerian radiation and radiation-hydrodynamic algorithms (Stone & Norman 1992a,b; Stone, Mihalas, & 
Norman 1992; Turner & Stone 2001; Hayes & Norman 2003; Hayes et al. 2005). This approach relies 
on staggered-mesh schemes to treat both the hydrodynamic and radiation components of the flow. In the 
original paper, Stone & Norman (1992a) lay out a hydrodynamics algorithm, which we extend to treat dense- 
matter hydrodynamics. More recently, Turner & Stone (2001) have extended this work to treat radiation- 
hydrodynamics in a comoving-frame, gray, flux-limited diffusion (FLD) approximation. The algorithm we 
describe in this paper extends this work to a multispecies, nonlinear, multigroup approach that is capable 
of treating non-LTE continuum problems. Ideally, we would like to solve the comoving-frame Boltzmann 
equation in multiple dimensions. However, the computational cost of such an effort would be so great 
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that multidimensional supernovae simulations could be carried out only for short times, even when using 
the most advanced, present generation of parallel-computing architectures. The multispecies, multigroup 
flux-limited diffusion (MGFLD) approach we detail here is an important stepping-stone to this goal. 

Finite-difference hydrodynamic algorithms possess several advantages that make them especially suited 
to core-collapse supernova simulations. First, these algorithms are formulated in a generalized, orthogo- 
nal coordinate system that allows for the easy interchange of Cartesian-, cylindrical-, or spherical-polar- 
coordinate systems. Second, these algorithms avoids use of Riemann solvers, either exact or approximate. 
This allows these algorithms to be extended to incorporate an arbitrary complex, and possibly non-convex, 
EOS. Additionally, such finite-difference methods can be easily extended to add new physics. Finally, these 
methods are relatively straightforward to implement on massively parallel, distributed-memory computing 
architectures. 

Our scheme similarly builds on this staggered-mesh approach, but makes major extensions to the treat- 
ment of radiation. This permits incorporation of numerous important aspects of neutrino flow. The first 
of these extensions accommodates multiple species of radiation (v e , v e , V M , V M , v T , v T ), in which particle- 
antiparticle types are coupled through pair-production processes. The particle-antiparticle coupling man- 
dates the simultaneous solution of the transport equations for both particles and antiparticles. The second 
extension is our multigroup treatment of the radiation spectrum, which includes ability to treat full coupling 
between energy groups that occurs by processes such as neutrino-electron scattering. Opacities and emis- 
sivities are included in an energy-dependent form, obviating need for any mean-opacity approximations. 
The third extension of capabilities is the treatment of Pauli-blocking effects in the radiation microphysics. 
This extension adds nonlinearities to the neutrino transport equations and leads to additional steps in the nu- 
merical solution of those equations. Our algorithm also allows modeling of matter deleptonization, through 
the addition of an equation describing evolution of electrons. Finally, our treatment of the hydrodynamics 
equations permits a complex EOS by incorporating a set of nonlinear solution algorithms for the Lagrangean 
portion of the gas energy equation. 

Much recent attention has been focused on the need for systematic processes of verification and vali- 
dation of computer simulation codes (Roache 1998; Knupp & Salari 2002; Calder et al. 2002, 2004; Post & 
Votta 2005). These procedures are required to achieve a reasonable degree of quality assurance of simula- 
tion results. Verification has been loosely defined (Knupp & Salari 2002) as testing to ensure that equations 
are being correctly solved, while validation has been defined (Roache 1998) as testing to ensure that mi- 
crophysical models are adequate descriptions of nature. In this paper, we present results of a number of 
verification problems that stress important components of our algorithm. Here, we do not concern our- 
selves with problem-specific microphysics. Validation tests for core-collapse supernova problems will be 
addressed in a forthcoming publication (Swesty & Myra, in preparation). 

An important consideration in the design of neutrino-radiation-hydrodynamics algorithms is that they 
be implementable on state-of-the-art massively-parallel computing architectures. We have designed our al- 
gorithm with this goal in mind, and have realized a parallel implementation in the form of a code (V2D), 
which we currently employ to simulate supernova convection (Myra & Swesty, in preparation) and proto- 
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neutron star instabilities (Swesty & Myra, in preparation). Although the focus of this paper is on the algo- 
rithm, not the implementation, we have aimed to provide all detail necessary to allow other developers to 
implement the algorithm and reproduce our results. 

This remainder of this paper is organized as follows: In §2, we introduce the coupled equations of 
neutrino-radiation hydrodynamics. Section 3 contains a description, in schematic form, of our algorithm 
for solving these equations numerically. (In appendices, we provide a detailed description of the finite 
differencing, numerical solution, and implementation of boundary conditions.) We present in §4 the results 
of verification tests that we have used to benchmark this algorithm. Finally, in §5, we present our conclusions 
about this algorithm. 



2. Equations of Neutrino Radiation Hydrodynamics 

The equations of neutrino radiation hydrodynamics must describe the time evolution of two primary 
components: matter and neutrino radiation. The matter is assumed at all times to be in local thermodynamic 
equilibrium (LTE). The neutrinos are never assumed to be in LTE, although such a situation may obtain in 
certain situations. For the moment, we assume that the radiation can be of an arbitrary form (e.g., photons, 
neutrinos, etc.), and we will make the distinction specific as needed. This allows our algorithm to be used 
for a variety of radiation-hydrodynamic situations that may, or may not, involve neutrinos. However, we will 
assume that multiple species of radiation are present. In the case where the radiation component consists of 
neutrinos, it is necessary to describe six different species of neutrino: v e , v e ,v^, V M , V T and v T . 

For the purposes of this paper, we assume the spatial domain to be free of macroscopic electric and 
magnetic fields. In principal, there is no reason that the algorithm we present here could not be extended to 
encompass magneto-hydrodynamic phenomena, but such extensions are beyond the scope of this work. 

In the subsections that follow, we first consider the hydrodynamic equations that describe the flow of 
dense matter, the equations that describe the evolution of the comoving multigroup radiation energy density 
in the flux-limited diffusion approximation, and the microphysical coupling between matter and radiation. 



2.1. Hydrodynamics with Neutrino-Radiation Coupling 

The starting point for a description of the material evolution is the set of Euler equations, which describe 
the dynamics of the matter. The corresponding starting point for a description of the radiation is the set of 
multigroup flux-limited diffusion equations, which we address in the next section. Since matter and radiation 
do not evolve independently, there are coupling terms that appear in both sets of equations to describe the 
transfer of energy, lepton number, and momentum between matter and radiation. The Euler equations, for 
the system under consideration, are: 

^+V-(pv)=0, (1) 
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%+V-(n e v)=N, 



(2) 



— + V-(£"v)+PV-v + Q: Vv = S, 



(3) 



d(pv) 



+ V-(pvv)+VP + V-Q + pV<£ + V-P rad = P. 



(4) 



dt 



Equation (1) is the continuity equation for mass, where p is the mass density, and v is the matter velocity. 
These quantities, and those in the following equations, are understood to be functions of position x and time 
t. Equation (2) expresses the evolution of electronic number density, where n e is the net number density 
of electrons over positrons. It is only relevant to include this equation if there is a variation in the ratio of 
the number density of electrons to the number density of baryons within the spatial domain, or if there are 
processes that can change the net number of electrons in the system. Thus, if the radiation being considered 
is electromagnetic, equation (2) is a redundant linear multiple of equation (1). However, in the presence of 
weak interactions in dense matter, equation (2) is usually independent of equation (1), and its right-hand side 
is non-zero. Here, we express that right-hand-side term — the net number production rate of electrons, having 
dimensions of number per unit volume per unit time — by N. To conserve lepton number, such reactions also 
imply a net number production of radiation from neutrinos or some other lepton. Therefore, evaluation of N 
involves integration of production rates over all neutrino energies and a summation over all neutrino flavors 
for any electron-number changing weak reactions (see §§2.3 and 2.4). The detailed microphysics of such 
reactions have been explored elsewhere, including Fuller et al. (1982); Fuller (1982); Fuller et al. (1985); 
Bruenn (1985); Langanke et al. (2003); Hix et al. (2003, 2005) and are beyond the scope of this paper. 

Evolution of the internal energy of the matter is given by the gas-energy equation (3), where E is 
the matter internal energy density, P is the matter pressure, and Q is the viscous stress tensor. Again, the 
right-hand side of this equation is non-zero whenever energy, of any sort, is transferred between matter 
and radiation. This can occur with neutrinos as a result of weak interactions or with photons as a result 
of electromagnetic interactions. For the moment we lump all such exchanges into the quantity S, which 
has dimensions of energy per unit volume per unit time, and represents the net transfer rate of energy 
from radiation to matter. The reactions that comprise S depend on the physical phenomena being modeled. 
However a general form for these reactions that encompasses most situations will be delineated in a later 
section of this paper. In the case of photons a detailed description of such reactions can be found in Mihalas 
& Mihalas (1984); Castor (2004); Pomraning (2005). For the case of neutrinos, in addition to the references 
for neutrino number-changing reactions listed above, additional reactions have been studied by Beaudet, 
Petrosian, & Salpeter (1967); Yueh & Buchler (1976, 1977); Schinder & Shapiro (1982); Bruenn (1985); 
Schinder et al. (1987); Mezzacappa & Bruenn (1993c); Ratkovic et al. (2003); Dutta et al. (2004), among 
others. 

Finally, equation (4) is the gas-momentum equation, where <1> is the gravitational potential, P ra d is 
the radiation-pressure tensor, and P is the net transfer rate of momentum due to microphysical interactions 
between radiation and matter. 



In equations (3) and (4), we have followed Stone & Norman (1992a) in our addition of a viscous 



-6- 



dissipation tensor to the Euler equations in order to account for dissipation that occurs in shocks. The details 
of of this viscous dissipation tensor are discussed in Appendix E. 

We note that it is also possible to substitute for equation (3) linear combination of the gas energy and 
gas momentum equations to get an evolution equation for the total matter energy (cf. Mihalas & Mihalas 
1984, eq. 24.5). 



For core-collapse supernova simulations, however, an internal energy formulation, as given in equation (3), 
has advantages. This is because there is a vast amount of internal energy in matter relative to kinetic energy. 
This follows from the thermodynamic domination of degenerate electrons, which contribute a large amount 
of zero-temperature energy and pressure. Given this situation, our choice of solving the gas-energy equation 
helps insure an accurate calculation of entropy, which is critical in degenerate regimes where a small change 
in energy can lead to a large change in temperature. In other problems, such as high mach-number flows, 
where kinetic energy dominates, a system may be better solved by using equation (5). 

Closure of this system of equations requires additional relations. First, is an equation of state (EOS), 
which is a parametric description of gas pressure and internal energy in terms of temperature, density, and 
composition. We discuss this further in §3.7. Second, is an expression for the gravitational potential <!>, 
which is discussed in Appendix F. Finally, microphysical expressions are needed to evaluate N, §, and P. A 
general form for these terms is discussed in §2.4 and Appendix I. 



For simulations of neutrino radiation-hydrodynamic phenomena, solutions of the full discrete ordinates 
Boltzmann equation, including comoving-frame and group-to-group coupling terms, have been attempted 
only in one spatial dimension (Mezzacappa & Bruenn 1993a,b,c; Liebendorfer et al. 2001, 2004). This is 
because of the computational cost associated with the high dimensionality of the Boltzmann equation. For 
simulations in more than one spatial dimension, the computational burden of solving the Boltzmann equation 
currently necessitates resort to an approximate solution. The recent 2-D work of Livne et al. (2004) ignored 
coupling terms to achieve computational simplicity and the 2-D work of Buras et al. (2006) used a variable 
Eddington factor method at low angular resolution to render the calculation tractable. 

In contrast, we implement a fully two-dimensional, Eulerian, multigroup, flux-limited diffusion scheme 
that keeps all order v/c coupling terms and which is practical for high spatial resolution simulations on 
current parallel architectures. This scheme extends our earlier work (Myra et al. 1987; Swesty, Smolarski, & 
Saylor 2004) as well as that of Turner & Stone (2001) and involves the solution of the zeroth angular moment 
of the Boltzmann equation. These equations take the form of a pair of angle-integrated, monochromatic, 
radiation energy equations in the co-moving frame that describe radiation of a particle and, where applicable, 




(5) 



2.2. Radiation Transport 
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its antiparticle: 

^ + V- (E e y) + V- F e - e-^ (P e : Vv) = S e , (6) 

-^-£+V-(£ e v)+V-F e -£— (P e :Vv) =S e . (7) 

These expressions are equivalent to equation (6.49) in Castor (2004), as derived by Buchler (1983). The 
scalar quantities E e and E e are the particle and antiparticle monochromatic radiation-energy densities at 
position x and time t. The particle and antiparticle monochromatic radiation-energy flux densities are given 
by vectors F e and F e . The particle and antiparticle monochromatic radiation pressure are given by P e and 
P e , which take the form of second-rank tensors. For the definition of all of these quantities we refer the 
reader to the comprehensive work of Mihalas & Mihalas (1984). The right-hand side quantities, § e and S e , 
account for coupling between matter and radiation. They contribute to the quantities N, § and P of equations 
(2)-(4). The form of this contribution is described in §2.4 and Appendix I. Expressions of the form P e : Vv, 
indicate contraction in both indices of the second-rank tensors P e and Vv. For photons, and other particles 
that are their own antiparticles, the barred expressions have no meaning and equation (7) can be ignored. 

Equations (6) and (7) actually represent a large set of equations. There is a pair of such equations 
for each wavelength or frequency in the spectrum of radiation. Additionally, if one is transporting more 
than one species of radiation particle {e.g., neutrinos of different flavors or some other collection of diverse 
particles), there will be additional sets of equations of this form to account for these additional species. 

Although the set of moment equations represented by equations (6) and (7) is exact, it does not possess 
a unique solution because of the multiple unknowns (radiation energy density, flux density, and pressure) 
in each equation. (However, note that in the hydrostatic limit, the terms involving the pressure tensor 
vanish.) The solution of the monochromatic radiation energy equation requires the specification of a closure 
relationship relating E e , F e and P e . Unless one has already solved the full Boltzmann equation (obviating the 
present discussion), the true relationships among these quantities are only known in the asymptotic limits 
of transport behavior — diffusion, where the optical depth is large; and free-streaming, where the optical 
depth is small. Therefore, solution of the monochromatic energy equation in general situations requires an 
approximate closure relationship. 

One of the most common approximations invokes the assumption that radiation is diffusive and obeys 
Fick's Law 

F e = -D e VE e . (8) 
In the diffusive limit, this Fick's law approximation becomes exact and the diffusion coefficient is given by 

where jcJ is the total opacity, given by 

jcJ = k£ + K° e + J deV(e,e'). (10) 
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This total opacity consists of contributions from the absorption opacity Kg, the total conservative scatter- 
ing opacity fc|, and the total non-conservative scattering opacity, which is expressed by the integral in 
equation (10). We use the subscript e to indicate that the quantity is a function of radiation wavelength, 
frequency, or energy. We discuss the absorption and non-conservative scattering opacities further in §2.3. 

However, the flux can become acausal if both Fick's Law and the diffusion coefficient of equation (9) 
are employed unconditionally. This is because F g is proportional to the VE e , which is unbounded. This 
problem usually manifests itself in optically translucent or optically thin (free-streaming) situations. Main- 
taining causality demands that |F e | < cE e always. One standard technique developed to maintain causality is 
flux-limiting, which has an extensive history associated with it (Minerbo 1978; Pomraning 1981; Levermore 
& Pomraning 1981; Bowers & Wilson 1982; Lund 1983; Levermore 1984; Cernohorsky, van den Horn, & 
Cooperstein 1989; Cernohorsky, & van den Horn 1990; Janka 1991, 1992; Janka et al. 1992; Cernohorsky 
& Bludman 1994; Smit, van den Horn, & Bludman 2000). In this paper, we will follow Myra et al. (1987) 
and Turner & Stone (2001) and make use of the flux-limiting scheme derived by Levermore & Pomraning 
(1981). However, our algorithm can easily be used, with virtually no modification, with other flux-limiting 
schemes that are based on the Knudsen number. 

In the Levermore-Pomraning closure, the flux is written in the form of Fick's law, but modifications are 
made to the diffusion coefficient to insure that causality is maintained and correct physical behavior occurs 
in the free streaming limit. A general form of a flux-limited diffusion coefficient is given as 

D « = -^> (ID 
which becomes the Levermore-Pomraning specification by defining the flux-limiter k e (R e ) as 

X e {R e ) = 2+Re (12) 

v ; 6 + 3R e +R e 2 

The quantity, R e , is the radiation Knudsen number, which is the dimensionless ratio of the radiation mean 
free path to a representative length scale. Thus, the Knudsen number is given by 

S -Wi <13) 

where the e subscripts emphasize that all these quantities are dependent on the energy of the radiation under 
consideration. This definition of the Knudsen number ensures correct limiting behavior, both when R e — > 
in the diffusive limit and R e —> oo in the free streaming limit. 

Irrespective of any approximation, the tensorial Eddington factor, X g , which relates radiation pressure 
and energy, is defined by 

P e =X e £ e . (14) 
The quantity X g is often written in the form of another general expression, 



X e = i(l-2 e )l + i(3* e -l)nn, 



(15) 
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where I is the identity tensor and where nn is a dyad constructed from n, the unit vector parallel to the 
radiative flux. The quantity % e is referred to as the scalar Eddington factor. Upon applying the Levermore- 
Pomraning prescription, we obtain the useful expression, 

% e =X e {R e ) + {X e (R e )} 2 R 2 , (16) 

which gives us the full Eddington tensor, 

X e = l - (l - X E {R e ) - {X £ {R e )f R e 2 ) I + \ (3Xe(R e ) + 3{X e (R e )} 2 R e 2 -l) nn, (17) 

in the Levermore-Pomraning scheme. 

With the application of flux-limited diffusion, closure relations are now entirely determined, and all 
moments of radiation are expressed in terms of E e . In similar fashion, a closure scheme for equation (7) 
follows by direct analogy. Thus, equations (6) and (7) can be cast in a form in which they possesses, at least 
formally, a unique solution, 

£^{(X e £ e ):Vv} = S e , (18) 

£ ^{(X e £ e ):Vv} = S e . (19) 
It is this form of the transport equation for which we describe a solution method in §3. 



dt 

d£ e 
dt 



+ V-(£ e v)-V-(D e V£ e )- 



+ V-(£ e v)-V-(D e V£ e ) 



2.3. Collision Integral 

The right-hand side of equation (6), the collision integral, can be expressed in a general particle- and 
species-independent way as 

= Nemis-abs + Npairs + Pdscat ■ (20) 

These terms account for various mechanisms by which energy may transfer between matter and radiation. 
Most radiation processes fall into one of these three forms and can be included in our algorithm. 

The first term on the right-hand side of equation (20), [Se] em i s -abs' represents emission-absorption of 
radiation by processes that change the monochromatic radiation energy or number densities. In photon 
transport, an example of such a process is the transition of an atom between different energy states that 
results in emission of a photon. In neutrino transport, an example is the capture of an electron by a nucleus 
that results in emission of a neutrino. In general, these processes can be expressed as 

Nemis-abs = ^ (* + H^e) " CK$E e , (21) 

where S e is the emissivity of the radiation field (with dimensions of energy per unit volume per unit time per 
radiation-energy interval). It is the rate at which energy is added to the radiation, while K" is the absorption 
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opacity for the reverse process (in units of inverse length). By making the flux-limited diffusion approxima- 
tion we have made the assumption that the distribution function in the collision integral is isotropic and thus 
the expression aE e /e 3 is the quantum-mechanical phase-space occupation number for the radiation field at 
position x, time t, and energy e. The quantity a is given by (he) 3 /4ng = 9.4523 MeV 4 cm 3 erg -1 for both 
photons and neutrinos. This follows from the statistical weight factor, g, being unity for both particles. 

The factor 77 takes on different values, depending on the quantum statistics of the radiation field under 
consideration. It is unity for photons and all other bosons, leading to the well-known stimulated emission of 
photons. For neutrinos and all other fermions, r\ = — 1, leading to inhibited emission — a term we find more 
physically intuitive than stimulated absorption (Bludman 1977), which is frequently used in the literature. 
This form follows naturally from the Pauli exclusion principle, which allows only a single fermion per 
quantum state. In the case of neutrinos, once the Fermi sea is fully occupied the emissivity drops to zero. 
Bruenn (1985) gives a more complete description of the quantum mechanical origin of this factor in the case 
of neutrinos. Finally, for a classical radiation field, r\ = 0, reflecting the Maxwell-Boltzmann character of 
classical particles. 

The quantum mechanical principle of detailed balance requires that S e and fc| be related. When ra- 
diation and matter are in chemical equilibrium, we have must have a relationship between emission and 
absorption (the forward and inverse reactions) such that the right-hand side of equation (21) vanishes. This 
balance relationship is expressed in Kirchoff 's Law, 



When radiation is in chemical equilibrium with matter, it makes sense to assign the radiation a chemical 
potential, which we represent by jx e . A chemical potential obviously has little meaning in non-equilibrium 
conditions; however, Kirchoff's Law always has meaning for determining the microphysical relationship 
between emission and absorption under any radiative conditions, whenever the matter is in LTE (Mihalas 
& Mihalas (1984), p. 387). To satisfy the detailed balance requirement when matter and radiation are out 
of equilibrium, one substitutes for jj, e in equation (22) the value it would take if matter and radiation were 
already equilibrated. This allows Kirchoff's Law to set the relationship between S e and fcf . The correctness 
of this procedure is a consequence of detailed balance, which must be satisfied microscopically, regardless of 
any macroscopic state of the system. As an example of its application, in the weak charged-current reaction 
e _ +^->v e +n, we substitute ;ii e = \i e + \l p — \l n , where H e ,p,n w& the electron, proton, and neutron chemical 
potentials for the matter in LTE. Note that if we were to apply this procedure for a process emitting photons, 
jU £ = at all times. 

The second term on the right-hand side of equation (20), [Se] pa i rs , represents processes that create a 
particle-antiparticle pair. (Since a photon is its own antiparticle, this process also represents production of a 
pair of photons.) In photon transport, an example of such a process is the mutual annihilation of a positron- 
electron pair to produce a pair of gamma rays. A corresponding analogy in neutrino processes is e~ + e + 



S e =cB e K a e {l-r]e^-^l T ), 



(22) 



where B e is the generalized Planck "black-body" function given by 




(23) 
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annihilation to produce a neutrino-antineutrino pair. In general, these processes may be expressed in the 
form 

Npairs = + (l + n^E E ) e J de'G(e,e>) (l +1^') , (24) 

where G(e, e') is the pair-production kernel for production of a particle with energy e and antiparticle energy 
e'. The monochromatic energy density of the produced antiparticle is given by E e i. In fermion transport, the 
r\ = — 1 factor gives two final-state Pauli blocking terms, reflecting the inability of a pair-production process 
to produce a fermion-antifermion pair in which either of the presumed final states is already occupied. In 
contrast, for bosons, final-state degeneracy is actually enhanced, since tj = 1. 

We note that equation (24) makes no provision for an inverse reaction. Although an inverse radiation 
pair-annihilation reaction is of potential importance under some conditions, and its addition to the algorithm 
is straightforward, we do not consider it in this paper. Radiation pair-annihilation reaction rates are usually 
strongly dependent on the angular distribution of the radiation and this situation is fundamentally incompat- 
ible with the assumption of isotropy that was made in deriving the monochromatic radiation energy equation 
from the Boltzmann equation. Inclusion of pair-annihilation terms thus requires some ad hoc assumption 
about the angular distribution of radiation that is dependent on the particular phenomenon being modeled. 
Detailed discussion of the role of pair annihilation will appear in future work on core-collapse supernovae 
where such effects are potentially relevant. 

The final term on the right-hand side of equation (20), [§ e ] scat , represents general, non-conservative, 
scattering processes. These processes result in no net creation or destruction of particles. Hence, their 
contribution to N is zero (see §2.4). However, they will change the energy distribution of the radiation field 
as a result of interactions with matter. 

For non-conservative scattering processes, the contribution to the collision integral is 

Nscat = ( 1 + T 7^) C / rfeV ( e > e ')£e'-£ e c / deV(e'.e) (l+n^E^ , (25) 

where k s (e,e') is the scattering opacity for particles in energy state e scattering into energy state e'. In 
analogy to the other processes discussed in this section, there is also final-state enhancement (or blocking) 
in these expressions for bosons (or fermions) resulting from the 1 + rjaE e /e 3 terms. 

Finally, we note that for each interaction presented in this section, there is possibly a conjugate reaction 
of importance that entails interactions involving antiparticles. These produce [S e ] versions of each of the 
terms we have presented. To calculate these interactions, one substitutes "barred" versions for each of the 
production and opacity terms, and interchanges E e and E e . In analogy with equation (20) the sum of all such 
contributions yields S e . The sole exception to this is equation (24) where the antiparticle analog is given by 

PMp- = + + V* 8 ) £ I d£ ' G ^'^ { l+T] ^ E£ ) ■ (26) 

Note that the same pair-production kernel G appears in equations (24) and (26) with the order of the argu- 
ments reversed between the two equations. 
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2.4. Radiation-Matter Coupling 

A key step in closing equations (l)-(4) is providing a method for evaluating the right-hand-side col- 
lision terms that couple matter and radiation. The evaluation of the sources for (N, §, and P) in terms of 
the results of §2.3 is straightforward. The right-hand side of equation (2) gives the rate of lepton exchange 
between the radiation an matter. It can be written as 

N = -£|i('S e -'S e )je, (27) 

where the form of the emissivity, '§ e , is given by equation (20) and its specifics by the subsequent equations 
in §2.3 for each flavor of radiation. The leading superscript / is used to denote the flavor of the radiation, 
e.g., for neutrinos / = e, }i, or T. The integral over e accounts for contributions from the complete spectrum 
of the radiation field. The minus sign in equation (27) accounts for the fact that a gain in lepton number 
for the radiation field is a loss for lepton number in matter. The sum over I accounts for the possibility of 
multiple species of radiation that can engage in number exchange with the matter. In practice, the number 
exchange is always due to electron neutrino-antineutrino emission-absorption. The barred term is non-zero 
if there are distinct antiparticles being evolved separately from the particles (e.g., both electron neutrinos 
and antineutrinos). Note that scattering is not a number-exchanging interaction and, hence, the contribution 
of the scattering terms, when integrated and summed, is zero. 

Regardless of whether there is a number quantity that is exchanged between matter and radiation, there 
is generally a non-zero energy exchange between the two. The net result of this exchange on the matter side 
is given by, S, the right-hand side of equation (3), which can be written as 

§ = -£|( f S e +%) de. (28) 

It is important to note that unlike equation (27), the sum is now over i (rather than /), which is a summation 
over all species of radiation, not just those that exchange net number with the matter. Once again, the barred 
term is non-zero whenever antiparticles are being evolved distinctly. 

Finally, we will ignore momentum exchange, P, between matter and radiation in our present flux- 
limited diffusion scheme, i.e., 

P = 0. (29) 

The calculation of microphysical momentum exchange between matter and radiation requires a knowledge 
of the angular distribution of the radiation or, at least, the knowledge of the angular averaged absorption 
opacity. In the case of neutrino transport in core-collapse supernovae, 1-D Boltzmann simulations have 
shown this effect to be negligible, and thus we ignore it for the remainder of this paper. However, it would 
be easy to include this effect for photons or neutrinos if the angular distribution of radiation were known. 
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2.5. Enforcing the Pauli Exclusion Principal for Neutrinos 

The equations of neutrino radiation hydrodynamics described in the previous subsections are semi- 
classical in that the Pauli exclusion principal for neutrinos is taken into account only in the collision integral 
terms. The left-hand-side of equations (18) and (19) are purely classical and do not guarantee that the 
occupancy of a specific neutrino energy state is less than or equal to unity. This constraint can be stated 
mathematically as 

< f e = ^E e < 1 (30) 

where f e is the neutrino distribution function and where where a = (he) 3 /Ang = 9.4523 MeV 4 cm 3 erg" 1 
for both photons and neutrinos (for which g = 1). The numerical solution of equations (18) and (19) for 
the neutrino and antineutrino spectral radiation energy densities can produce values of E e and E e for which 
the distribution functions have values greater than unity. This is most likely to occur in situations where the 
neutrino distribution function becomes highly degenerate, such as in the core of a proto-neutron star. This 
problem has been known for some time (Bruenn 1985). 

Since distribution function values greater than unity are obviously unphysical for neutrinos we need 
to supplement the solution of the neutrino radiation-hydrodynamic equations by adopting an "enforcement" 
algorithm that is ensures that the constraint represented by equation (30) is satisfied after a new value of the 
neutrino spectral energy densities is computed. We detail this enforcement algorithm in Appendix L. 

2.6. Conservation 

The hydrodynamic equations that we have presented in this section are not in a conservative form that 
would admit a finite-volume approach to their discretization. Therefore, the equations do not guarantee 
exact numerical conservation of either energy or momentum. No discretization scheme can simultaneously 
conserve all physically conserved quantities. An excellent example of such quantities are linear and angular 
momentum, which cannot, in general, both be conserved by the same discretization. The neutrino radiation- 
diffusion equations we have presented are also not in conservative form. Some attempts have been made to 
arrive at a discretization of neutrino transport equations (Liebendorfer et al. 2004), in the 1-D Boltzmann 
case, that conserves both neutrino energy and number. However, in general no such discretization has been 
discovered. In order to ensure that sufficient accuracy is being achieved in a simulation, one must monitor 
the conservation of various physical quantities that are important. The relative importance of conservation 
is problem dependent and we will address this for this algorithm in the context of core-collapse supernovae 
in a future work. 



3. The Numerical Method 

The algorithm that we employ for the solution of the 2-D monochromatic radiation-hydrodynamic 
(RHD) equations, described in the previous section, is an extension of the ZEUS family of algorithms of 
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Stone & Norman (1992a,b); Stone, Mihalas, & Norman (1992); and Turner & Stone (2001). The extensions 
include the incorporation of a multigroup treatment of the radiation spectrum with group-to-group coupling, 
the incorporation of multiple radiation species with particle-antiparticle coupling and Pauli blocking, the 
solution of an additional advection equation describing the evolution of the electron number density, and 
the incorporation of a complex equation of state. We also employ a methodology for solving the implicitly 
discretized radiation-diffusion equations that is different from the approach set forth by Turner & Stone 
(2001). 

The algorithm we employ solves the hydrodynamic equations (1), (2), (3), and (4), along with flux- 
limited diffusive transport equations (18) and (19). Our algorithm is Eulerian and employs a staggered 
mesh similar to other ZEUS-type algorithms (Stone & Norman 1992a,b; Stone, Mihalas, & Norman 1992; 
Turner & Stone 2001; Hayes & Norman 2003; Hayes et al. 2005). Like these other algorithms, our time 
evolution scheme utilizes a combination explicit techniques to evolve the hydrodynamic portions of the 
RHD equations while employing implicit techniques to solve the transport portions of these equations. Our 
time evolution scheme is dissimilar to these other schemes in that the order of solution of substeps differs 
from these schemes. 

In the subsections that follow, we describe the generalized computational mesh that we employ, the use 
of operator splitting, the order of operator-split substeps, the discretization of the equations, and the parallel 
implementation of the method. 

3.1. Computational Mesh 

We employ a spatially staggered mesh on an orthogonal coordinate system identical to that employed 
by the ZEUS family of algorithms (Stone & Norman 1992a,b; Stone, Mihalas, & Norman 1992; Turner & 
Stone 2001; Hayes & Norman 2003; Hayes et al. 2005). For the sake of simplicity, we have not allowed the 
mesh to adapt dynamically in any way and we consider the mesh as fixed in time. It would be straightforward 
to extend the algorithm we describe in this paper to accommodate the moving mesh that is described in Stone 
& Norman (1992a). 

We number cell edges in each of our two coordinate directions, x\ and X2, by integer indices. The ith 
cell edge in the x\ direction has an x\ coordinate given by [jci],-, while the j'th cell edge in the X2 direction 
has an X2 coordinate given by [x2]j. The cell centers have coordinates in each direction given by [jci],-+(i/2) 
and [jC2]y+(i/2)- Thus, the location of a cell-centers in this mesh is fully specified by a pair of discretized 
coordinates in thexi and ^2 directions ([jti]j+(i/2)> [ x 2\j+(i/2))- This staggered mesh is illustrated in Figure 1. 

Intensive quantities, such as pressure, mass density, internal energy density, temperature, electron num- 
ber density (or electron fraction), etc. are defined at cell centers. Components of vector quantities, such ve- 
locities, momenta, gradients of intensive quantities, and fluxes are defined on the corresponding cell faces. 
We refer to these latter quantities as face-centered variables. The spatial location for each of these types 
of variables is depicted in Figure 1. Note that we use standard subscript notation to define the discretized 
analogs of all quantities, e.g., [T] i+ ^/ 2 )j+(3/2) is the discretized temperature variable defined at coordinates 
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([ x i]i+(i/2)> [ x 2];+(3/2))- Our finite difference notation is described in Appendix A. 

It is clear from Figure 1 that vector components are not co-located at a single spatial point on the 
mesh. Occasionally, quantities derived from these components are needed at alternate locations, in which 
case a scheme that averages values from nearby spatial locations is employed to compute values at the 
needed location. We discuss such averaging on a case-by-case basis when we detail the discretization of the 
equations. 

In order to discretize the spectral variables, we also define a mesh over the energy dimension, i.e., 
the spectrum of radiation energies. The range of the domain in the energy dimension is discretized into 
groups, i.e., cells in energy space. The kth group has a lower edge with an energy coordinate [e]k and the 
center of the kth group has an energy coordinate given by [s]k+(i/2)- Discretized spectral quantities, such 
as spectral radiation energy densities and spectral flux densities, are usually defined at group centers. Since 
such quantities usually have a dependence on spatial location as well as energy, the discretized analogs of 
these quantities carry an extra subscript indicating their location in the energy dimension. For example, the 
radiation energy density E e at spatial coordinates ([jci];+(i/2)j [^2]y+(i/2)) an d energy coordinate [e]fc+(i/2) 
is denoted as [£'e]yt+(i/2).(+(i/2),y+(i/2)- The location of energy-dependent intensive quantities such as the 
spectral radiation energy density, and energy-dependent vector quantities, such as the spectral flux density 
are illustrated in Figure 2. A complete listing, delineating where various physical quantities are defined in 
the spatial-energy meshes, is found in Appendix B. 

3.2. Covariant Formulation 

We choose to follow Stone & Norman (1992a) by writing all finite-difference expressions in terms of a 
generalized orthogonal coordinate system that is capable of describing Cartesian, cylindrical, and spherical- 
polar coordinate systems. Our goal is to enable a single code that is easily adaptable to any 2-D curvilinear 
coordinate system, avoiding the labor that would otherwise be required to implement a code in each individ- 
ual coordinate system desired. This technique is well described in (Stone & Norman 1992a) and we refer 
the reader there for details. The notation for coordinates and other geometrical quantities that we employ 
in each coordinate system are described in Table 4, located in Appendix C. The detailed form of the metric 
coefficients, the gradient and divergence operators, and tensor expressions that are needed to evaluate the 
radiation-hydrodynamic equations are described in their entirety in Appendices H and J. 

3.3. Operator Splitting 

Our algorithm employs operator splitting to decouple the overall time integration of the radiation- 
hydrodynamic equations into substeps. The motivation for this procedure is discussed in Stone & Norman 
(1992a), to which we refer the reader. In general, we split the right-hand-sides of the time evolution equa- 
tions into advective, source, viscous, and radiation-matter-coupling terms and solve these split equations to 
update the hydrodynamic and radiation quantities accordingly. 
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[ x 2]/+(l/2) 




[ x 2]j-(l/2) 



Fig. 1 . — A portion of the staggered spatial mesh, showing the location of spatial coordinates. The bold 
black lines define cell edges while the gray lines show the coordinates of cell-centers. The location of a 
typical cell-centered quantity, y, and of a typical vector quantity, G, with components Oi and 02, are also 
shown. 
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Fig. 2. — The full three-dimensional mesh in x\-X2-£, showing the positions for evaluation of the spectral 
radiation quantities within a cell. These quantities are evaluated at points on the plane that passes through 
the mid-point of the grid in the radiation-energy dimension. The face-centered radiation variables {e.g., 
spectral radiation flux density) are evaluated at the positions shown by the black or green dots on the cell 
faces. The cell-centered variables {e.g., the spectral radiation energy density), are evaluated at the cell center, 
the position of which is shown by the red dot. The blue dots, which lie on cell faces in the energy direction 
show the position of fluxes representing transfer of energy between energy groups. 
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The following describes the application of this operator-splitting approach to the equations in our 
model. The time integration of the continuity equation (1) requires no operator splitting, since there is only 
a single advective term, and no source or collision term, in the equation. Thus, we can restate equation (1) 

as 
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The electron conservation equation (2) is split into two terms, 
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is the source or collision-integral term. In a similar manner, the gas-energy equation (3) is split into four 
separate sets of terms: advection terms, the Lagrangean or source terms, viscous dissipation terms, and the 
collision-integral terms. 
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The gas-momentum equation (4) is operator split into five sets of terms 
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where the advection term is 



the source terms are 



the viscous dissipation terms are 



the radiation pressure terms are 



and the collision integral terms are 
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Finally, the radiation-energy equation (18) is operator split as 
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while the combination of the diffusive and collision integral terms are defined by 
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The antiparticle monochromatic diffusion equation (19) is operator split in the analogous fashion to equa- 
tion (18). 

We also note that each of the aforementioned advection terms is itself directionally operator split into 
two pieces corresponding to advection in each of the two coordinate directions which we generically denote 
as x\ and X2- Thus for each set of advection terms we can write 
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dt 
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+ 



dt 



(50) 



Because of the complexity of the operator split equations, we restrict the discussion of the numerical 
methods used to solve the individual pieces of the operator equations to Appendices D-L. In the next 
section we concentrate on the order of updates, based on these operator split pieces, employed to evolve the 
equations from time [t] n to time [t] n+l . 
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3.4. Order of Solution of Operator Split Equations 

Our algorithm employs the following order for solution of the operator-split equations detailed in the 
previous subsection. The complete sequence of solving each of these operator split pieces constitutes the 
algorithm for evolving the equations from time [t] n to time [t] n+l = [t] n +At. A schematic illustration of 
our algorithm, for a single timestep, is provided in Figure 3. The details of each substep are provided in 
Appendices D-L. The hydrodynamic portions of the operator-split equations are solved explicitly, while the 
radiative portions of the equations are solved implicitly The motivation for this choice of a hybrid explicit- 
plus-implicit approach is well-described in Stone, Mihalas, & Norman (1992), and we refer the reader there 
for a detailed description of the issues involved. 

Following Stone & Norman (1992a), we denote partial updates of variables of quantities by means of 
superscripts. Thus, at the beginning of a timestep, the matter internal energy density is denoted by [E] n . The 
partially updated internal energy resulting from, for example, substep / (updated via eq. [40]) is denoted 
by [E] n+ f . The superscript n + f serves to indicate that the quantity includes all partial updates prior to and 
including substep /. The final update of each quantity, within a given timestep, is labeled by superscript 
n+\. When we denote a discretized quantity without spatial- or energy-index subscripts, we are referring 
to the entire spatial or energy range of that discretized quantity. In our description of each substep, we will 
describe what quantities are updated as a result of that substep. 

The substep a-b in the algorithm consists of explicit numerical solution of the advection portions of 
the operator-split equations. This substep is actually a combination of two directionally split substeps cor- 
responding to advection in the x\ and X2 directions. In this substep, all advective portions of the operator 
split equations are solved, namely, equations (32), (34), (37), (42), and (48). Note that, since the radiation 
energy density is a function of spectral energy, equation (48) must be solved for every energy group and for 
each type of radiation present, i.e., for each of the six types of neutrinos. Thus, equation (48) represents a 
set of 6N g equations that must be solved, while equations (32), (34), (37), and (42) represent five equations 
(the momentum equation is actually two equations, one for each of the the two components of the veloc- 
ity). Thus, when the number of energy groups N g is large, the computational cost of the advection substep 
is dominated by the cost of solving the radiation-advection equations represented by equation (48). Our 
algorithm for explicit solution of the advection portion of the operator-split equations is exactly the same 
as that of Stone & Norman (1992a) and is detailed in Appendix D. We note that after the calculation of 
the updated values of the radiation energy densities, the Pauli exclusion principal constraint enforcement 
is applied. This is indicated in Figure 3 by red shading of the advective update box. We also note that 
the advective update substep is itself composed of numerous substeps in which the advection equations are 
solved by directionally split substeps. The directionally split advection algorithm utilizes Norman's consis- 
tent advection scheme (Norman et al. 1980) in which the advection of all quantities is tied to the mass-flux. 
The details of the advection substeps are illustrated in Figure 28 (see Appendix D), to which the reader is 
referred for more detail. The net result of this substep are the partially updated quantities [p]" +/ \ [E] n+h 
[T} n+b , [P] n+b , [n e ] n+b , [Y e ] n+b , [\} n+b ,[s] n+b , [ e E e ] n+b , [ e E e ] n+b , [»E e ] n+b , [»E e ] n+b , [ r E e ] n+b and [ T E e ] n+b . 

The second and third substeps (substeps c and d) in the solution of the operator-split equations involve 
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Fig. 3. — The algorithm for advancing the model by a single timestep from [t] n to [t] n+l . The red boxes 
indicate steps where the Pauli exclusion principal constraint is enforced after new values of the neutrino 
energy densities are calculated. The variables listed in each box are those that result from the update (those 
in parenthesis are inputs to the update step). "Enter" signifies the beginning of the timestep, while "exit" 
signifies the end of the timestep. 
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the evolution of equations describing the radiative evolution of electron neutrinos and antineutrinos and the 
exchange of energy and lepton number between matter and these neutrinos. This substep involves the im- 
plicit solution of the set of radiation diffusion equations for electron neutrinos and antineutrinos represented 
by equation (49) and the collision-integral equations represented by equations (35) and (40). In the second 
step, represented by the second box (c) in Figure 3, the complete set of implicitly differenced diffusion 
equations for electron neutrinos and antineutrinos represented by equation (49) is solved simultaneously via 
Newton-Krylov iteration. The details of the finite-differencing and numerical solution of these equations 
is detailed in Appendix H. Once the implicit solution of this set of equations has been accomplished, the 
amount of lepton and energy exchange between matter and electron neutrinos and antineutrinos is fixed. Af- 
ter the new values of the electron-neutrino-antineutrino energy densities are calculated, the Pauli exclusion 
principal constraint enforcement algorithm is applied to the electron-neutrino-antineutrino energy densities. 
The application of this constraint-enforcement algorithm is indicated by the red shading of the second box 
(c) in Figure 3. The second substep results in updated quantities [E] n+C , [T] n+C , [P] n+C , and fully updated 
radiation quantities [ e E e ] n+l and [ e E e ] n+l . 

In the substep d (represented by the third box [d] of Figure 3), since the amount of energy and lepton 
exchange with matter has been fixed by the previous substep, equation (35) is solved for the new value 
of electron number density and, thus, the new value of electron fraction Y e . Subsequently equation (40) 
is solved implicitly for the new value of internal energy density. Once the new internal energy density is 
determined, the equation of state determines the new value of matter temperature and pressure corresponding 
to the new internal energy density. The details of this substep are described in Appendix H. The third substep 
results in the fully updated quantities [E] n+d , [T] n+d , [P] n+d , [n e ] n+d , and [Y e ] n+d . 

The second and third substeps (substeps c and d) are subsequently repeated for the muon neutrinos 
and antineutrinos in substeps e and / (shown as boxes e and / the flowchart) and tauon neutrinos and 
antineutrinos in substeps g and h (shown as boxes g and h in the flowchart). In substeps e and g, the Pauli 
exclusion principal constraint algorithm is applied to the muon neutrino and antineutrino energy densities 
and the tauon neutrino and antineutrino energy densities, respectively. This is indicated by the red shading of 
the boxes corresponding to substeps e and g in Figure 3. In substeps / and h, the solution of equation (35) is 
not required since the production of muon neutrinos and antineutrinos and tauon neutrinos and antineutrinos 
results in no change in lepton number — these neutrinos are always produced in particle-antiparticle pairs. 
Equation (40) is solved for the new matter internal energy density, temperature, and pressure, as described 
in the case of the second substep. Substep e results in the updated quantities [£ , ]" +e , [r]" +e , [P]" +e , [^E e ] n+ , 
and [PE e ] n+l . Substep / results in the updated quantities [E] n+ f, [r] n+ ^, and [P] n+ f. Substep g results in the 
updated quantities [£]"+* [T] B+ *, [P] n+g , [ c E e ] n+ \ and [ c E e ] n+1 . Substep h results in the updated quantities 
[E] n+h , [T] n+h , and [P] n+h . 

In substep i, the momentum and velocities are updated via the solution of equation (43) to account for 
gravitational- and pressure-induced accelerations. This substep is almost identical in detail to that of Stone 
& Norman (1992a), but we describe this in detail in Appendix F. In this paper, we consider the gravitational 
force to be spherically symmetric based on the mass constained interior to a given radius. The description of 
the calculation of the gravitational mass is also detailed in Appendix F. This substep results in the updated 
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quantities [v] n+i and [s]' l+! '. 

In substep j, the momentum and velocities are updated via the solution of equation (45) to account 
for radiation-pressure-induced accelerations. This substep relies on the Eddington factor differencing of 
the gray transport algorithms of Stone, Mihalas, & Norman (1992) and Turner & Stone (2001) which, for 
our multigroup transport, is applied on a group-by-group basis. We described this approach in detail in 
Appendix J. This substep results in the updated quantities [v]" +; and [s]" +; . 

In substep k, the components of the artificial viscous stress are calculated according to the prescription 
of Stone & Norman (1992a). This calculation is described in Appendix E. 

In substep £ the momentum and velocities are updated via the solution of equation (44) to account for 
accelerations induced by the gradients of the viscous stress. This substep is identical in detail to that of 
Stone & Norman (1992a), and we describe this in detail in Appendix E. This substep results in the updated 
quantities [v] n+1 and [s] n+1 . 

In substep m, the internal energy density is updated via the solution of equation (39) to account for 
viscous dissipation. Like substep I, this substep identical in detail to that of Stone & Norman (1992a), and 
we describe the update equation in Appendix E. This substep results in the updated quantity [E] n+m . We 
point out that the temperature and pressure are not updated in this step, as they will be updated after the 
following step. 

In substep n, the Lagrangean portion of the gas energy equation, described by equation (38) is solved to 
account for compression or expansion of the gas and the effects of viscous stresses. The time differencing of 
this equation is implicit. However, the since the divergence of the velocity in equation (38) is evaluated based 
on the partially updated velocities [\] n+l , there is no spatial coupling between the unknowns in equation (38). 
This equation can thus be solved by a local, nonlinear iterative solution algorithm in each spatial zone. The 
finite differencing of equation (38) and our solution algorithm are described in Appendix G. This substep 
results in the updated quantities [Ef + \ [Tf + \ and [P] n+l . 

This sequence of partial updates represented by substeps a, b torn constitutes the algorithm for evolving 
the equations of neutrino radiation hydrodynamics from time [t] n to time [t] n+l = [t] n + At. 



3.5. Boundary Conditions 

Up to this point, we have neglected any discussion of boundary conditions and how they are applied 
within the algorithm. In general, boundary conditions for a specific quantity are applied immediately after 
any update of that quantity. Thus, any given quantity may have boundary-condition updates several times 
during the course of a single timestep. The details of how specific boundary conditions are applied are 
delineated in Appendix K, and we refer the reader there for more information. 

In a parallel implementation of our algorithm, where parallelism is achieved via spatial domain decom- 
position, there are also internal "processor boundaries." Values of variables at these boundaries must be kept 
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consistent among processors. Thus, boundary updates are a frequent occurrence, since internal processor 
boundaries are in the middle of actively computed regions. Such consistency requires update of values of 
a quantity at the edge of processor boundaries after each update of that quantity — before it is needed for 
another calculation requiring spatial derivatives of that quantity. We discussion this issue in §3.8. 



3.6. Timestep Selection 

The stability of our algorithm is restricted by the stability of the solution of explicitly solved operator- 
split equations. In calculating the timestep we follow the algorithm laid out in Stone & Norman (1992a). 
This algorithm for selecting the timestep depends on several different types of stability criteria which are 
then combined as an RMS average to yield a stable timestep. 

The calculation of the timestep is based on four key criteria. The first is the Courant timescale (calcu- 
lated in both the x\ and x 2 directions), which represents the minimum sound-crossing time for a particular 
zone in each dimension. More formally, 

[A?Co U rant]-+( 1 i/2) J+ (l/2) = T min ([^+(1/2) > M j+(l/2)) > ( 51 ) 

where c s is the local speed of sound at coordinate ([jci]/+(i/2), [ x i\j+{\/2))- An accurate calculation of c s 
requires the equation of state to supply an adiabatic index. For the purposes of timestep selection, however, 
a conservative overestimate of c s is obtained by making a conservative approximation {i.e., overestimate) of 
the sound speed by using a polytropic EOS having an overestimate of y. 

The second and third metrics are the flow timescales in the x\ and X2 directions, which are the timescales 
over which a particle in the fluid, located at one cell face, can travel to the opposite face, in each respective 
direction. These timescales are expressed as 

I- *i fl°wJ«+(l/2)j+(l/2) - mm 7^1 ' TT] ( 52 ) 

and 

r A/ ,»+! _ • ( [g2] ;+(1 / 2 ) [Ax 2 ] J+(1/2) [g 2 ] i+{1/2) [Ax 2 ] j+{l/2) \ 

[At X2 flowJi +( i/2)j + (i/2) - mm nr] » n~] • ( 53 ) 

Finally, the fourth timescale is a viscous dissipation timescale set by the magnitude of the viscous stress. 
This timescale is defined as 

lAt r+ i _ 2 rr^J ^W/2) [ g2 ]jA, 2 ] y+(1/2) \ 

^convJ,- + (i/ 2 ),y + (l/2) - Z V l q mm 777] FT71 ' T7T] rTH ' W 

Vi U lJ/+lJ+(l/2)-[ U lJiJ+(l/2) lV2\i +{ l/2),j+l-[V2\ i+{ l/2)J J 

where we compare the timescales in each direction and take the minimum. The quantity l q is a number of 
order unity and is defined in Appendix E. 
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The inverse squares of each of these timescales are added at each mesh point. The maximum value of 
the quantity in parenthesis in equation (54) is found for the entire spatial domain. The inverse square root 
of this quantity represents the minimum timescale in the entire domain. A fraction, which we refer to as the 
CFL fraction and designate by /cfl, of this time is used as the new timestep size. It is our practice to fix 
/cfl throughout the course of a given simulation. Typically, we set / C fl = 0-5- 

Thus, the timestep value eventually used is 

I 1 1 

max j 1 2 

all 'J y [A?Courant] ; -+(l/2)J+(l/2) faxi flow ]; + (i/ 2 )j+( 1/2) 



[A?* 2 flow] ; - + (i/2)j+(i/2) [^conv];+(i/ 2 )j+(i/2) 





3.7. Equation of State and Opacity Interface 

This algorithm makes no assumption about the equation of state other than assuming that the EOS is 
of the form 

E = E(T,p,Y e ) (56) 

and 

P = P(T,p,Y e ). (57) 

Our numerical algorithm can accommodate an arbitrary EOS of this form. The algorithm does not rely on 
solution of a Riemann problem and makes no assumptions about convexity in the EOS. Thus, the EOS can 
be supplied as a simple formula that can be evaluated in a small subroutine or, alternatively, in the more 
general form of a thermodynamically consistent tabular interpolation (Swesty 1996). The algorithm does 
require that it be possible for the relationship described by equation (56) be inverted, either analytically or 
numerically to yield 

T = T(E,p,Y e ). (58) 

In particular, whenever a new value of the matter internal energy E is calculated, it must be possible to 
compute a new value of the temperature T corresponding to that internal energy. A polytropic EOS can be 
easily accommodated within the relationships described by equations (56)-(58). 

We also assume that the absorption opacity, the conservative scattering opacity, and the emissivity be 
of the form 

K a e = K a e (T,p,Y e ,e), (59) 
4 = 4(r,p,y e ,£), (60) 

and 

S e =S e (T,p,Y e ,e). (61) 



-26- 



The absorption opacity and the emissivity should be related in such a manner as to preserve the quantum 
mechanical principle of detailed balance (see §2.3). 

The scattering opacity is assumed to be of the form 

K-l e , = Kl el (T,p,Y e ,e,e'), (62) 

where, in the scattering reaction, e is the energy of the incoming particle and e' is the energy of the outgoing 
particle. This opacity should also preserve detailed balance. 

Finally, the pair-production rate is assumed to be of the form 

G e , e , = G e , e ,(r,p,y e ,e,e'), (63) 
where e is the energy of the neutrino that is produced and e' is the energy of the antineutrino that is produced. 



3.8. Parallel Implementation 

The size of problem encountered in multidimensional radiation-hydrodynamic models, particularly 
in stellar collapse, necessitates our use of massively parallel computing resources in order to solve the 
discretized equations. Since we solve a long-timescale problem, it is necessary that we achieve strong- 
scaling, i.e., we wish a fixed-size problem to scale well to a large number of processors and, therefore, 
reduce our wall-clock time to solution. In addition, the number of variables in the problem requires a large 
amount of memory, further necessitating a parallel solution strategy. 

A parallel implementation of our algorithm can be achieved via spatial domain decomposition of the 
2-D spatial domain into a logically Cartesian topology of 2-D subdomains. Such a decomposition is illus- 
trated in Figure 4. This is a well-established scheme for the parallelization of numerical methods for partial 
differential equations (Gropp et al. 1999a). The logically Cartesian process topology is straightforward to 
create using MPI (Message Passing Interface) subroutine calls (Gropp et al. 1999b) and can be configured 
to allow for periodic boundaries if so desired. This logically Cartesian spatial decomposition is independent 
of the choice of coordinate system and is carried out with an orthogonal spatial mesh defined by generalized 
coordinates, which we have described previously. The partitioning of this mesh into subdomains is accom- 
plished by specifying the size of the process topology in each of its two dimensions. Once the size of the 
process topology has been specified, the mesh is divided in an approximately even fashion over the process 
topology to achieve a balancing of computational work. By specifying the process topology to be as square 
as possible, the ratio of the number of ghost zones to non-ghost zones can be minimized, thus reducing the 
communication-to-computation ratio and improving the scalability of the algorithm. 

In our parallelization scheme, we do not decompose in the spectral energy dimension of our mesh. Thus 
the 2-D quadrilateral subdomains illustrated in Figure 4 are actually 3-D hexahedra, where the third dimen- 
sion is the spectral energy dimension. By not decomposing the problem in the energy dimension we avoid 
costly communication that would be associated with the evaluation of the integral terms in equations (24), 
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Fig. 4. — A schematic of the outer edges and center of the entire computational mesh, showing spatial 
domain decomposition into a logical Cartesian mesh for parallel computing. Blue shaded zones represent 
genuine boundary zones of the computational mesh. Zones that are shaded red represent internal ghost zones 
that are duplicate copies of zones that are resident on other processors. The actual size of the subdomains 
is problem dependent. In the present implementation of the algorithm the width of the boundary and the 
ghost zones is always two. This is a result of the present differencing scheme, in which coupling between 
cells within a given equation never extends beyond next-nearest neighbor. The x\ and %2 starting and ending 
indices on each subdomain are labeled as i s , i e and j s , j e respectively. 
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(25), (27), and (28), as well as the application of the Pauli exclusion principal constraint enforcement algo- 
rithm. 

Under a logically Cartesian spatial domain decomposition, our discretization algorithms for the hydro- 
dynamic and radiation-transport equations require a limited set of communication patterns. Evaluation of 
fluxes or spatial derivatives gives rise to a local discretization stencil that requires information from "ghost 
zones" surrounding each subdomain. The values of variables in these ghost zones must be obtained from 
adjacent subdomains by means of message-passing to and from nearest-neighbor processes. This process of 
ghost-zone exchange is illustrated in Figure 5. Ghost zone values of a specific quantity, such as density, must 
be exchanged before those values are needed in the evaluation of any expression in which those variables 
appear. These exchanges can be accomplished asynchronously, but we avoid discussion of the complexities 
of doing this, since that topic is beyond the scope of this paper. 

The number of the ghost zones required is a function of the discretization scheme. An examination of 
the finite difference equations, presented in the appendices, indicates that the maximum number of zones 
that couple within a single equation is five in each spatial dimension. For example, equation (1), for van Leer 
advection of scalars, couples five zones i— |, i — ^, i + 5, i + 1, and i+ | centered about the i+j cell center 
in the x\ direction. If a spatially higher-order differencing scheme were implemented, there would be longer- 
range coupling among zones. Therefore, the width of the ghost-zone region would be correspondingly larger. 
Readers who desire a more complete description of ghost-zone exchange are referred to Chapter 4 of Gropp 
et al. (1999a). For the implementation of the algorithm described in this paper, two ghost zones in each 
spatial dimension are required. 

Unfortunately, nearest-neighbor message passing is not the only communication pattern required by the 
algorithm. Global reduction operations are required in two instances. First, calculation of the timestep size, 
At (see §3.6), requires a global reduction to determine the minimum CFL time for the entire domain. Sec- 
ond, the Newton-Krylov-based iterative solution of the implicitly differenced radiation-transport equations 
(Appendix H), requires global reduction operations to evaluate vector inner products. In the BiCGSTAB 
algorithm, which we use to solve the linear systems in the Newton-Krylov scheme, this can mean multiple 
global reductions per iteration. This can impose a bottleneck to scalability for simulations running on large 
numbers of processors. To reduce this bottleneck, we have developed an algebraically equivalent variant of 
the BiCGSTAB algorithm, which requires only one global reduction per iteration (Swesty 2006, in prepara- 
tion). This improvement can be seen in Figure 6, where we plot the parallel speedup of the algorithm when 
calculating a supernova model on seaborg, the IBM-SP at NERSC. The major floating point cost of the 
Newton-Krylov algorithm is expended in the evaluation of the finite-differenced nonlinear diffusion equa- 
tion. This operation requires only nearest-spatial-neighbor communication to evaluate the finite-difference 
stencil of the divergence operator. Whenever the nearest neighbor is a zone whose data is resident on another 
processor, we amortize the communication cost by performing the nearest neighbor-communication asyn- 
chronously. This allows us to carry out, simultaneously, the portions of the matrix- vector multiply operation 
corresponding to interior (local) zones of each subdomain. This yields a further improvement in scalability, 
as seen in Figure 6. 



Fig. 5. — A schematic showing how logically adjacent processes can exchange internal boundary infor- 
mation through the exchange of ghost zones. Ghost zones contain data from the computational domain 
that is resident on a processor other than the one being considered. The mutual exchange of this internal 
"boundary" information ensures that each process in the calculation is working with up-to-date information. 
The arrows show the flow of information. Active computational zones in one processor's domain (e.g., the 
red-shaded zones in the bottom subdomain) are transferred to the ghost zones in the other processor's sub- 
domain. The green-shaded zones show the corresponding flow in the opposite direction. In the figure, the 
blue-shaded zones represent ghost or genuine boundary zones that are not involved in the exchange. The 
flow of information to and from these zones is governed by exchange rules for the transverse dimension. 
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Fig. 6. — Parallel speedup of the algorithm, with fixed problem size, on the NERSC IBM-SP (seaborg). 
Note that the roughly 82 Gflop/s, which is obtained running with 512 processors, represents at speedup of 
about 14 over a 32-processor run — a parallel efficiency of about 85%. 
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3.9. Code Structure 

Although time advancement is the most important part of the algorithm, it is only the central portion of 
a sequence of many operations that manage an overall simulation. For completeness, we present a diagram 
of the remaining portion of the algorithm in Figure 7. Several of the initial operations, depicted by blue boxes 
in Figure 7, reflect a parallel implementation of the algorithm. The main computational effort described by 
the flowchart is the time evolution loop depicted by the red boxes in Figure 7. The core of this loop is the 
timestep evolution previously described in Figure 3. The purpose of the loop in Figure 7 is simply to evolve 
the equations forward in time. 

Of note in Figure 7 is that we periodically write out data from a run in the form of checkpoint files. 
These files capture the state of the simulation and serve two functions. They act as restart files, so that a sim- 
ulation can be resumed from any checkpointed timestep. They also serve as a resource for post-processing 
analysis. Use of parallel file systems, MPI-I/O (Gropp et al. 1999b), and HDF5 2 in our simulations allows 
us to create a checkpoint file, using a particular number of processors and processor topology, and restart it 
at a later time, using a different processor count or topology. In fact, it is routine for us to use checkpoint 
files in this manner. The file portability built into both MPI and HDF5 also allows our use of these files 
across diverse computer architectures. 

4. Code Verification Tests 

To test our algorithm, we have subjected the code that implements it, V2D, to a suite of verification 
test problems that stress individual components of the code in a variety of different contexts. These tests 
are broken out into five main classes: (1) hydrodynamic tests, which test only hydrodynamic portions of 
the code, i.e., without radiation transport; (2) gravity tests, which also involve hydrodynamics, but which 
stress self-gravity of the system; (3) transport tests, which test the radiation transport portion of the code, 
but in a static medium (i.e., without hydrodynamics); (4) radiation hydrodynamics tests, which test the 
coupled hydrodynamic and radiation portions of the code; and (5) correctness and scalability tests, which 
are a diverse set of tests used to test such things as the nonlinear solvers implemented within the code and 
the correctness and scalability of the parallel implementation of the code. 

Because of the length of this manuscript, we have chosen only a few problems for each of these cate- 
gories. These problems were chosen to encompass the most important aspects of algorithm correctness. In 
some categories, such as hydrodynamics, there are many classic verification test problems that we have not 
included for reasons of brevity. All of the test problems we have included in this section, along with many 
others, are run as automated regression tests that are executed daily against our source code repository. 

In the following subsections we provide verification test descriptions and results for each of the afore- 
mentioned categories. 



2 http://hdf.ncsa.uiuc.edu/HDF5/ 
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Fig. 7. — An overview of the complete algorithm from program initiation to finish, focusing on the steps 
required to set up the problem. The timestep advancement algorithm, which is the red-colored loop, is 
detailed in Figure 3. 
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4.1. Hydrodynamic Tests 

The hydrodynamic verification tests we have implemented employ solely the hydrodynamics aspects 
of the algorithm. These problems only involve equations (l)-(4), as presented in §2.1, with all radiation- 
coupling terms omitted. 

We perform three standard tests of hydrodynamics. The first is a passive advection test, which tests the 
ability of an advection scheme to preserve features in a the profile of a physical quantity when the quantity 
is advected by a moving medium. The second and third problems, the shock tube and the Sedov-Taylor blast 
wave, test a code's ability to reproduce the solutions to two standard problems whose solutions are known 
analytically. 

Two of these tests, the passive advection problem and the shock tube problem were performed by Stone 
& Norman (1992a), but since our hydrodynamic algorithm differs slightly in the order of solution steps, we 
repeat these tests to establish the performance of our algorithm on the classic problems. The Sedov-Taylor 
blastwave problem was not addressed by Stone & Norman (1992a), and so we include the problem here. 

4.1.1. Passive Advection 

The passive advection test exclusively exercises advection equations described in Appendix D. The 
goal of this test is to delineate how faithfully a waveform is preserved as it propagates through a moving 
medium. Hence, it serves as an excellent measure of the diffusivity of the advection algorithm. For com- 
parison purposes, our setup closely follows the passive advection test results reported in Stone & Norman 
(1992a). Source terms in the hydrodynamics are ignored and only the advection equations are evolved in 
time. There are no gravitational forces, no pressure gradients, and no radiation. In Cartesian coordinates, 
the test uses a uniform, time-independent velocity field, which is set to 5 x 10 4 cm s _1 , parallel to one of 
the axes. As an initial waveform, we set up a square pulse 50 zones wide near the left-hand boundary of a 
domain 400 zones wide, where each zone is 1 cm wide. This pulse is applied to the density, and since we as- 
sume the material is an ideal gas, it also appears in the pressure and internal energy. Specifically, for passive 
advection in the x\ direction, we set p = 1 g cm~ 3 and P = 1 erg cm~ 3 in the region 5 cm < x\ < 55 cm, and 
p = 10~ 10 g cm~ 3 and P = 10~ 10 erg cm~ 3 elsewhere. The velocity in the x\ direction is set to 5 x 10 4 cm s _1 
and the velocity in the X2 direction is set to zero. The matter internal energy density is set to E = P/(y— 1) 
everywhere with 7=5/3. These initial conditions are evolved until the pulse propagates a distance of five 
times its initial width. With the background velocity, and using a timestep of 5 ;lls, this is achieved in 1000 
timesteps. 

Results for our test are shown in Figure 8. Although there is some diffusion evident, given that the 
square pulse shape is not preserved exactly, the performance of V2D's implementation of advection is 
in good agreement with the results of Stone & Norman (1992a), who have also implemented van Leer's 
scheme. We perform numerous variations on this test. In one such Cartesian variation, we run the mirror 
image of the test in the negative x\ direction. After making the necessary changes to adjust for the symmetry, 
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Fig. 8. — Results in Cartesian coordinates for passive advection of a 50-zone wide pulse carried along in 
the x direction for 1000 timesteps — a distance corresponding to 5 pulse widths. The advection scheme 
used is that of van Leer. We show comparisons of the following quantities: (a) mass density, (b) velocity, 
(c) pressure, and (d) internal energy density. Exact results, which would be expected from a "perfect" 
advection scheme, are shown in solid red, whose vertical lines are located at 255 and 305 cm. Numerical 
results are plotted with diamonds. In the velocity plot (b), the numerical data has been set to agree precisely 
with the exact velocity and, hence, completely obscures the analytic line. 
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results of the tests in the +x\ and —x\ directions are found to be bitwise identical. We also perform the same 
pair of tests in the +X2 and — %2 directions. Once again, results of this pair are bitwise identical, mutatis 
mutandis. This ^-direction pair also has bitwise identical results to the pair of tests in the x\ direction. 

4.1.2. Shock Tube 

The shock tube problem serves as an important test of both the advection scheme employed by a code 
and of the overall performance of a hydrodynamic algorithm. First introduced by Sod (1978), it is now 
widely used as a standard verification test. Although essentially a one-dimensional problem, it turns out 
to be a useful test for a 2-D code. This is because, in addition to checking for the correct behavior in the 
principal direction of interest, it can also check that this correct behavior is exactly replicated at all points in 
the second dimension. 

We set up and run the shock tube problem in each of the three coordinate systems (Cartesian, cylindri- 
cal, and spherical) and in each direction. For the purposes of this discussion, we choose ^ = x\ , where x\ 
can be related either to the x-direction in Cartesian coordinates, the radius 9? in cylindrical coordinates, or 
the radius r in spherical coordinates. We construct the numerical solution with 100 zones in the x\ direction 
and 12 zones in the X2 direction. As just discussed, the x\ coordinate is centered about zero. It extends 
4 cm in total such that — 2 cm < x\ < + 2 cm. Zones in each direction are uniform in spatial extent. 
The initial conditions are those set forth by Sod: in the left half of the domain we have p = 1 g cm~ 3 and 
P = 1 erg cm~ 3 , while in the right half we have p = 0.125 g cm~ 3 and P = 0.1 erg cm~ 3 . All velocities 
are initially zero. The internal energy density is initialized to E = P/{y— 1) everywhere. In this, and all the 
shock tube problems in our test suite, we have set the polytropic index, y = 5/3. 

In our tests, we follow the subsequent evolution until t = 0.70 s. The decision of how long to run the 
problem differs among authors. For comparison purposes, the particular choice is important only in that the 
problem needs to run long enough so that the resulting shock wave can propagate a significant distance into 
the low-density gas. However, the calculation must not run so long that the shock reaches the right-hand 
boundary of the computational domain. Our choice of timestep size is governed by the CFL condition (see 
§3.6), and our fraction of the CFL time is set to 0.5. With these choices, each test run takes 59 timesteps to 
reach 0.70 s. Results for a test run using spherical coordinates are shown in Figure 9. In all our tests, we 
use our implementation of the van Leer advection scheme combined with Norman's consistent advection, as 
described in Appendix D. Our result compares favorably to Stone & Norman (1992a) {cf. their Figure 11). 
The similarity of these results is not surprising since the hydrodynamic algorithms of ZEUS and V2D have 
a similar approach, and both calculations use van Leer advection. Our results also compare well to the best 
results in Hawley et al. (1984b) (cf. their Figures 6-15). In particular, the results using their consistent 
advection scheme are comparable to those shown here. 

In addition to running the shock tube problem in spherical coordinates, as shown in Figure 9, we also 
perform the shock tube problem in cylindrical and Cartesian coordinates. In the case of the latter, we perform 
the tests in the +x and — x directions and in the +y and — y directions. After adjustment for the different 
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Fig. 9. — Results for the 100-zone, non-relativistic, shock-tube problem in spherical coordinates, with the 
shock propagating outward in the radial direction. To assist in comparison of these results with those of 
Cartesian shock tubes, we define a radial coordinate x\ such that x\ = is the initial position of the shock 
wave. With this definition, x\ =r — 2 cm, where r is the distance from the center of the spherical "tube". We 
show comparisons for the following: (a) mass density, (b) velocity, (c) pressure, and (d) specific internal 
energy. These results compare favorably with those of Stone & Norman (1992a) and the best schemes in 
Hawley et al. (1984b). 
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coordinate systems and directions, the results of these tests are compared and found to be bitwise identical 
to each other. As a check of the parallelization of the hydrodynamics algorithm, we also run the shock tube 
problem in spherical coordinates in the +r direction on several different processor counts and topologies: a 
serial (lxl) single-processor run and parallel 4-, 16-, 64-, and 256-processor runs in 2x2, 4x4, 8x8, and 
16 x 16 topologies, respectively. When the results of these runs are compared, they are found to be bitwise 
identical to each other. 



4.1.3. Sedov-Taylor Blast Wave 

The Sedov-Taylor blast wave is another example of an idealized problem for which a self-similar an- 
alytic solution exists. The solution was originally found by Taylor (1950) and Sedov (1959). The problem 
is as follows: At time t = 0, a large amount of energy is deposited at a single point in an otherwise uniform 
cold, diffuse, ideal gas. The result is an expanding spherical shock front, behind which is an expanding, 
spherically symmetric, distribution of matter. Since the solution to the problem is self-similar, there exist 
similarity variables such that radial profiles of the physical quantities are invariant in time with respect to 
the similarity variables. 

For the numerical test, we set up the problem in spherical coordinates. We choose a spatial domain 
that is 1-cm in radius and is spanned by 100 uniform width zones in the x\ (radial) direction. In the X2 
(angular) direction, the domain size n radians with 6 equally spaced angular zones (the number of angular 
zones should be truly irrelevant to the outcome of the problem). A solution with multiple angular zones 
is performed purely as a check that the hydrodynamics code is capable of maintaining perfect spherical 
symmetry. The initial conditions are p = 0.1 g cm~ 3 and T = 10~ 8 MeV. We initialize the problem in two 
different ways. In the first case we place a high energy density "bomb" of T = 1 MeV, p = 0.1 g cm~ 3 
matter in the the innermost three zones. We refer to this case as the the three-zone bomb case. In the second 
case, the bomb is located only in the centermost zone and the temperature is set to T = 27 MeV so as to 
keep the total bomb energy constant. We refer to this as the one-zone bomb case. The precise values of the 
initial conditions are somewhat arbitrary as long as the vast majority of the energy in the problem is initially 
confined to a small region near the center. We assume an ideal gas with y = 5/3. We evolve the solution 
until a time of 7 x 10~ 9 seconds using a CFL factor of 0.5 which takes 1496 steps for the one-zone bomb 
case and 458 steps for the three-zone bomb case. The results for both cases are shown in Figures 11 and 
10. Our results for the three-zone bomb case agree well with other codes where the bomb is spread out 
over multiple zones. However, the one-zone bomb case is perhaps a more authentic way of initializing the 
problem which is rarely seen in the published literature. As we can see from Figure 11, the agreement of 
the numerical and analytic solutions is substantially worse in the one-zone bomb case. The one-zone bomb 
result can be improved somewhat by using non-uniform zoning. 

Defects caused by zoning in numerical solutions of the Sedov-Taylor problem are common. The prob- 
lem is especially difficult because most of the material piles up behind the shock, and it is difficult to 
maintain the spatial resolution at all the points necessary to track accurately its movement. Even codes that 
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Fig. 10. — Results for the three-zone bomb Sedov-Taylor blast wave in spherical coordinates plotted versus 
the self-similar analytic solution. All quantities are plotted in terms of the dimensionless similarity variables 
for the problem, as indicated in the legend. (The "j" subscript on the quantities in the ratios indicates that 
they are analytic values at the shock front.) The solid lines represent the analytic solution; the data points 
are the numerical solution. 



-39- 




Fig. 11.— 



Same as the previous figure except the bomb is confined to one zone. 
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can draw on adaptive mesh refinement (AMR) still have issues. For example, Fryxell et al. (2000), using the 
FLASH code, manage to track the position of the shock very well, but are less successful in obtaining good 
post-shock profiles for the quantities of interest. Pen (1998), also employing an AMR-based code, is forced 
to choose between optimizing either resolution and positioning of the shock or resolution of the post-shock 
material. 

Finally, we note that when we check our results for spherical symmetry, we observe that the solutions 
of all radial rays are bitwise identical to each other, and that all velocities in the angular fe) direction are 
identically zero at all times. 

4.2. Gravity Tests 

Although V2D is a code intended for calculations in two spatial dimensions, self gravity is currently 
approximated by assuming that the mass distribution interior to a point is spherically symmetric (see Ap- 
pendix F). This treatment naturally leads to gravity verification tests that have spherical geometry. The two 
gravity verification tests that we consider in this paper (a unit density sphere and a polytrope of index 1), 
both have spherical geometry and are solved in spherical coordinates. 

4.2.1. Unit-Density Sphere 

In this test problem we compute the gravitational potential for a sphere of unit density. 

Since the gravitational potential given by equation (F9) is a trivial function of the gravitational mass 
interior to a point, this problem reduces to a test of our algorithm for computing the mass interior to a 
point. In a serial implementation, this summation could be accomplished trivially. However, in a parallel 
implementation, this summation requires a series of reduction operations across the process topology. In 
practice, we accomplish this by first summing the mass in a shell at a given radius by a series of MPI global 
reduction operations in the X2 direction. Once the amount of mass at a given radius is known, the total mass 
interior to each radial zone is summed in an outward fashion starting from the inner edge of the grid. These 
summation steps involve numerous message-passing operations and the purpose of this test also serves to 
verify the correctness of our implementation. 

The problem is set up by constructing a 10-km sphere of unit density, which consists of 100 radial 
zones of equal spatial width and 40 angular zones of equal angular width, extending from 0-n radians. The 
gravitational mass is computed using the method given in Appendix F. The result is compared against the 
simple exact computation that is possible for this problem. 

The value of this test is the comparison that is possible between the exact method and the various modes 
of calculation that may take place when V2D is executed. One of the pitfalls in parallel computing is the 
potential effect in global summation operations caused by floating point arithmetic being neither associative 
nor distributive. Thus, the answer to a summation operation that is distributed over multiple processors will, 
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in general, produce a result that is dependent on the number of processors used. Furthermore, results can 
vary slightly from run to run, even on identical numbers of processors, since MPI (which is commonly used 
for message passing) offers no guarantees about the order of arrival of inter-processor messages. Because 
of this, the order of individual operations may vary from run to run, possibly producing variable results. 

To estimate the magnitude of these effects, we calculate the gravitational mass in a number of different 
ways and compare each result against the exact result. The deviations of these various results from the exact 
result are shown in Table 1. In what is labeled "MPI with deterministic reduction" in Table 1, we have 
replaced the MPI_ALLREDUCE operation in the standard MPI library with one of our own construction. 
Our version is designed to guarantee an identical order of all critical arithmetic operations, regardless of 
processor count. As the table shows, the deterministic results agree well with the exact result. As indicated in 
the table, these deterministic results also agree with each other bitwise, i.e., to the full numerical precision of 
the system on which they were run. This is to be expected from the deterministic design of these calculations. 
Although determinism is a desirable feature, it is obtained at the cost of serializing certain key time-intensive 
portions of the computation. Therefore, it is not a practical method to employ in long-running production 
calculations. 

The results in Table 1, labeled "standard MPI," use the ordinary, standard-compliant MPI libraries for 
all message passing. Such calculations proceed with no regard to determinism of results. (This is also 
the mode in which we do production computations using the gravity solver.) The results from this method 
show, curiously, a smaller deviation from the exact result than shown by the deterministic calculations. More 
significantly, this deviation differs with processor count, as expected. However, comparisons with the exact 
result are so good, that it is inconceivable that the overall results of a large-scale, long-running simulation 
could ever vary in an important way with processor count, solely because of the magnitude of variations 
exhibited here. Machine epsilon is of the same order (~ 10~ 16 ) as the deviations — for the architectures on 
which both these tests and our production runs are performed. 

On the basis of these tests, and of others that exercise other parallelized portions of the code (e.g., 
parallel solvers, the SN tests — see §§4.5 and 4.6), we can proceed with confidence that parallelization alone 
cannot produce incorrect results (due to mistakes in parallelization) or deceptive results (due to answers that 
depend on processor count). 

Figure 12 shows a comparison of the the exact result and that calculated by the deterministic method 
using a single processor. As already noted, the agreement is excellent. 

4.2.2. Stability of a Poly trope 

The polytrope stability verification test involves both gravity and hydrodynamics. This test is initialized 
by numerically constructing a static polytrope. We then measure how well the static solution is maintained 
when the polytrope is subjected to hydrodynamic evolution. This test provides a check on both the spherical- 
gravity solver and on the correctness and stability of the hydrodynamic portion of the V2D algorithm. 
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Fig. 12. — A test of V2D's spherical gravity solver. A 10-km sphere of unit density is constructed and 
the gravitational mass at each of 100 grid points is compared against an exact calculation. The results 
shown here are for a single-processor calculation corresponding to the first row of Table 1. The maximum 
discrepancy of the total mass calculated by the two methods produces relative differences of 2.4 x 10 . 
Plots of similar data, for the other versions of this test shown in Table 1 , are visually and (in the case of the 
deterministic runs) bitwise indistinguishable from this figure. 
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Table 1. Gravitational Mass Calculation Comparisons: Maximum Deviations of 
Various Test Runs vs. the Exact Result 



method used 


processor 


processor 


relative deviation 




count 


topology 3 


from exact method b 


MPI with deterministic reduction 


1 


lxl 


2.4 x 1(T 14 
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2x2 


2.4 x 1(T 14 




16 


4x4 


2.4 x 1(T 14 




64 


8x8 


2.4 x 1(T 14 




256 


16x16 


2.4 x 1(T 14 


standard MPI 


1 


lxl 


5.2 x 1(T 16 




4 


2x2 


4.1 x 1(T 16 




16 


4x4 


3.6 x 1(T 16 




64 


8x8 


2.6 x 1(T 16 




256 


16x16 


4.1 x 1(T 16 



a Processor topology dimensions refer to the logical structure by which the two- 
dimensional computational grid was decomposed among processors, not to any phys- 
ical processor-interconnect topology that the computer hardware employed. 

b The deviations shown here for each of the "MPI with deterministic reduction" 
results, in fact, agree with each other bitwise — they agree to the full numerical pre- 
cision of the system on which they were run. 



-44- 



We choose a poly trope with index n = 1, which corresponds to an polytropic equation of state that 
has an adiabatic index of 7 = 2. With this choice, we can construct a hydrostatic "star" by solving the 
Lane-Emden equation, which has an analytic solution for n = 1 (Chandrasekhar (1967), p. 84ff), 

p(^)=p c ^l, (64) 

where p c is the central density and t; is a radial coordinate defined such that the radius, 

r = /3£, (65) 

where 

The quantity K is the polytropic constant and G is the universal gravitational constant. As indicated by 
equation (66), for n = 1, /3 is independent of p c . For convenience, we set p c = 1 g cm~ 3 and choose 
K = 2nG, such that j3 = 1. Hence, for this problem, the density range throughout this configuration ranges 
between and 1 g cm~ 3 , and the surface is located at n centimeters from the center. 

This analytic solution yields initial conditions for the test problem. We construct a model having 100 
radial zones (equally spaced in radius) and 32 angular zones (equally spaced angularly). The zone center 
coordinates are determined by 

VUU2) = \yj {[Am? +[r]i + i[r]iH\r]if (67) 

in accordance with the coordinate regularization scheme of Monchmeyer & Miiller (1989). To initialize the 
density, we assign values for p at cell centers by averaging the analytic solution evaluated at the adjacent 
cell faces, i.e., 

1 / sin([*i],.) sin([xi] m )' 



[p] *-+(v^-+(V2) - 2 {-^t + ^tr) ■ (68) 

This initialization is carried out in an angularly symmetric fashion. Apart from the angular symmetry, we 
find the exact procedure by which the mesh is initialized to the analytic solution is unimportant to the 
subsequent evolution. 

Once initialized, the system is evolved through 1000 timesteps. The timestep size is chosen to have a 
CFL factor of 0.5, as determined by the procedure described in §3.6. After 1000 timesteps, this corresponds 
to a model time of about 14 500 s. (To prevent the necessity of tiny timestep sizes, given the relatively small 
angular zone size when using 32 angular zones, we do not include contributions to the Courant time in the 
angular direction from the first 10 radial zones. To ensure that unphysical motions do not occur because of 
this, we also suppress all angular motion, that is motion in the X2 direction, inside these 10 zones.) 

Figure 13 shows a comparison between the evolved configuration and the analytic solution. Except 
near the center, where grid resolution is an issue, the numerically evolved solution typically agrees with the 
analytic solution to a few parts in 10 4 or better. The larger deviation near the center is a result of gridding 
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Fig. 13. — The normalized density for the polytropic stability test after 1000 timesteps. Except near the 
center, where grid resolution is an issue, the numerically evolved solution typically agrees with the analytic 
solution to a few parts in 10 4 or better. 
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effects where the mesh provides insufficient spatial resolution. The radial velocities are also checked. After 
1000 timesteps, the absolute value of the radial component of velocity nowhere exceeds 9 x 10~ 7 cm s _1 . 
The transverse component of velocity remains identically zero at all timesteps throughout the spatial domain. 
The previous statement is the equivalent of saying that the code allows all radial rays to evolve in identical 
fashion — as expected, and as required. Similar results are obtained for the case where we choose K = 
3.874 x 10 4 cm 5 s~ 2 g _1 and the mass of the configuration is M = 1.4M . 



The transport verification tests are intended to establish the numerical performance of the radiation 
diffusion aspects of the algorithm. Each of the tests described in this section assume that radiation is propa- 
gating through a hydrostatic medium. Thus, the hydrodynamic variables are not evolved. Problems that test 
coupling of radiation diffusion to hydrodynamics are incorporated in §4.4. In all transport and radiation- 
hydrodynamic tests, the Pauli blocking factors that appear in the collision integral are retained with tj = — 1. 
However, the problems are such that the particle occupancies are much, much less than unity and therefore 
the blocking factors obtain a classical value of unity. 



This problem makes use of an analytic solution of the time-dependent diffusion equation in 1-D Carte- 
sian coordinates for the case of a constant diffusion coefficient (Crank 1980; Butkov 1968; Riley, Hobson, & 
Bence 1998). For initial conditions in the form a Gaussian profile centered at point ro, the analytic solution 
to the diffusion equation is given by 



where t Q is the initial time from which the problem is evolved, and 8t is the evolution time. The quantity 
E the initial amplitude of the pulse and D is the diffusion coefficient. This solution can be used to test 
the diffusion algorithm in both Cartesian and spherical coordinates. In the case of spherical coordinates the 
pulse is centered at some large radius such that the plane-parallel limit of the spherical coordinate system is 
reached. 

The initial conditions for the numerical solution are constructed as follows. We use a single energy 
group of radiation, bounded by and 1 MeV. We use the arithmetic average of the group boundaries, 
0.5 MeV, to define the energy at the group center. Since the problem is hydrostatic and there is no emission- 
absorption, there will be no change in spectral shape. Hence, the details concerning energy grouping are 
irrelevant. Spatially, we divide the domain into 100 radial zones, covering 1 cm radially, with each zone 
spaced uniformly. In order to achieve the plane-parallel limit the inner edge of the radial grid is placed at 
x\ = 10 4 cm. Flat (zero-flux) boundary conditions are applied at the inner and ±X2 edges of the grid with free 



4.3. Transport Tests 



4.3.1. Diffusion of a 1-D Gaussian Pulse 




(69) 
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streaming boundary conditions at the outer edge. The number of angular zones is arbitrary to the outcome 
of the problem (each radial ray should evolve virtually identically). We perform our test with 16 angular 
zones, equally spaced in angle, covering a range of 10~ 4 radians and check to make sure the solution in each 
radial ray is the same. This choice of angular coordinates makes the computational domain in the transverse 
direction cover about 1 cm and, therefore, gives the domain the shape of a small curvilinear "square." To 
allow a more complete testing of the code, we perform our calculation by setting up three identical parallel 
transport problems, thus testing all three species of radiation that V2D is currently set up to track. This 
is done for both particles and antiparticles of each species. In the language of neutrino transport, this is 
the equivalent of evolving electron, muon, and tauon neutrinos and their antiparticles. This makes I in the 
equations below have values from 1 to 6, inclusive. 

The discretized initial profile is given by 

I \r\. . /„> — r„ 

4Dt a 



[ £ eW(i/2),;+(i/2)j+(i/2) = ^exp <^ -— — ( [r] i+(1/2) - r ) } for all k,£,i,j, (70) 



where, as indicated, the radial profile is replicated over all angular zones, and applied to all species of 
radiation. Here, k = 1 for the single energy group of radiation being evolved. We have set the diffusion 
coefficient D = 10 7 cm 2 s _1 , which in the diffusive limit corresponds to a total opacity, Kj = 10 3 cm -1 via 
equation (9). The constant E Q = 1 erg cm~ 3 MeV -1 . After evolving for a time t, the analytic solution to this 
problem leads us to expect a numerical solution of the form 

1 &]* + (l/2),K-(l/2),y + (l/2) = E ^J^8 t ^{- 4 D{t + 8t) ( M i+(V2) " r °) } 

for all k, £,i, j, (71) 

where we have set t = 1.0 ns for this test. We run the simulation to where t + 8t = 4.3 ns, a total of 100 
timesteps. The timestep size was chosen to be At = 100 Ar/c. Figure 14 shows a comparison between 
the analytic solution and that achieved by the numerical test. Except near the center, where the numerical 
solution slightly lags the analytic solution, the agreement is excellent. With the setup described, and using 
the familiar Euclidean norm, the test exhibits a small relative residual at the final timestep, 

/ \2\ V 2 

V V 1 1 \J7 I analytic I \tt 1 numerical ' ' 

LiLj y L £ 'eJ / t + (l/2),,-+(l/2)J+(l/2) PeJ*+(l/2),;+(l/2)J+(l/2) 



' 2 = I * .w „ ^ , ,„.w , j _ am9) (72) 

analytic 

t+(l/2),i+(l/2)J+(l/2) 



ZiLj( e [Ee]T mC 



for each radiation species individually. As a check, we confirm that the results for each species are bitwise 
identical. 

We also check how well the code maintains spherical symmetry. Because the metric tensor, when us- 
ing spherical coordinates, contains trigonometric functions of the angular coordinates, it cannot be expected 
that spherical symmetry will be identically preserved. In our 16-angular-zone model, we find that the max- 
imum relative deviation in the radiation energy densities between corresponding points along radial rays is 
4 x 10~ 15 . Thus, spherical symmetry is well maintained, and the departure from it is of the degree expected 
for a system whose machine epsilon is of order 10~ 16 . 
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Fig. 14. — Results after 3.3 ns (100 timesteps) for a spreading Gaussian pulse of radiation, propagating in 
the radial direction in spherical polar coordinates. 
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4.3.2. Diffusion of a Gaussian Pulse in Two Dimensions 

This problem extends the 1-D Gaussian pulse problem to two dimensions and tests the performance 
of the diffusion algorithm in multiple dimensions. The computational domain is similar to the spherically 
symmetric problem: The same single radiation-energy group extending from 0-1 MeV is used. The spatial 
domain is a curvilinear square similar to that used in the previous problem, but this time measuring 2 cm 
on a side. The center of this domain is located at a radius of 10 4 cm from the origin in a spherical-polar 
coordinate system. The radial domain is divided into 101 zones; the angular domain of 2 x 10~ 4 rad is also 
divided into 101 zones. The use of an odd number of zones ensures that the pulse can be placed at the exact 
center of the domain. Since the computational domain is small relative to its distance from the origin, it 
resembles a Cartesian grid, which is the way we view it in Figure 15, which shows both initial conditions 
and results. Nevertheless, the numerical problem is actually solved in spherical polar coordinates and the 
solution is compared against the analytic solution in that same coordinate system. 

The initial Gaussian pulse is given by 

Ee=E ° eXp { Z ^blf L } (73) 

which gives an azimuthally symmetric distribution about the point r D , which is the displacement of the 
center of the pulse from the origin. For our test problem, we have selected r D = 10 4 cm and the direction of 
the displacement to be located at the angular center of our chosen domain. Once again, we have chosen a 
uniform, constant diffusion coefficient D = 10 7 cm 2 s -1 , E a = 1 erg cm~ 3 MeV -1 , and t = 1 ns. 

This pulse is initialized numerically at time t D to the values given by 

eW(l/2),/+(l/2)J+(l/2) = 
E ° eXP {~4bT (Wi+a/z) 008 ^ - [ e U+(i/2) " [ ];+(i/2)) - foCOsCF- [0] jt +(1/2) - e o )) } 

XeXP {"4zk (Mi+(i/2)Sin(^- [ U+(i/2) " Wj + (i/2))-roMV-[0]j s+ (i/2) ~ %)) 

for all k,£,i, j, (74) 

where *P is half of the angle subtended by the computational domain at the origin, [0]y s+ ( 1( / 2 ) is the angular 
coordinate at the first half-integer mesh point, and [0] 7 - + n m is the angular coordinate of a general cell center 
where the distribution is being initialized. The quantity O is the angular coordinate of the center of the initial 
pulse. Flat (zero-flux) boundary conditions are applied at the edges of the grid. We allow the initial pulse to 
evolve for 3.3 ns in time. As in the 1-D problem, the timestep was chosen to be At = 100 Ar/c. By direct 
substitution, it is straightforward to show (Riley, Hobson, & Bence 1998) that at time t + t, the analytic 
solution for this problem is 

£ '= £ °(^) MP {-4D(^7) |r - r ° |2 }- <75) 
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Fig. 15. — Initial conditions and solution for the diffusive propagation of a two-dimensional Gaussian pulse. 
For display purposes, we show everything in an approximate Cartesian coordinate system, though the actual 
calculations are carried out in spherical polar coordinates. The "Cartesian" coordinate x is defined such that 
x = r — r , where r = 10 4 cm. Similarly, the "Cartesian" coordinate y is defined as r Q. The initial Gaussian 
pulse is shown in (a), where both the colormap and the contours display the energy density over the spatial 
domain. In (b), the solution is shown after ~ 3.3 ns of evolution, using the same colormap as (a). Except for 
areas close to the boundary, where edge effects are present, the pulse precisely maintains its axisymmetric 
form about (x,y) = (0,0) throughout the evolution. 
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We note that for this two-dimensional problem, the fraction containing t and t, which multiplies the Gaus- 
sian function, appears as a linear power. This contrasts with the 1-D case, where this factor appears under a 
square root (eq. [71]). 

The analytic solution corresponds to a predicted numerical solution of 

1 [ £ eW(l/2),i+(l/2)j+(l/2) =E ° {j~^ t/ 

XeXP {~ 4D(l + t) (W^(l/2) COS ^- [ b' s + (l/2) - I W/2)) -r COs(»F- [0]y s+(1/2) " O ) 

i-^ ([r],. +(1/2) sin^ - [d] js+(l/2) - [0] J+m ) - r sm(V - [0] . +(1/2) - O )) ' 

for all k,£, (76) 



XeXP ^4D( f , 



Figure 15 shows the results of solving the test problem numerically with V2D. In Figure 15(a), we 
show the initial conditions as given by equation (74). Both the colormap and the plotted contours show the 
initial radiation energy density. 

Figure 15(b) shows results after ~ 3.3 ns of evolution using the same scaled colormap for E e as in 
(a). The initial pulse has spread. Except for areas close to the boundary, where edge effects are present, 
the pulse precisely maintains its axisymmetric form about (x,y) = (0,0) throughout the evolution. The 
boundary discrepancies result from our imposition of flat boundary conditions (see Appendix K) and our 
use (by necessity!) of a finite domain. This results in a numerical solution that is too small at the boundary. 
The behavior of the solution near the boundaries could likely be improved through the imposition of time- 
dependent boundary conditions. Nevertheless, at 0.1 cm inside a boundary, the deviation of the numerical 
from the analytic solution is already reduced to < 1.4%, which is typical of the deviation over the rest of 
the domain. The character of the deviation of the numerical from the analytic solution is similar to the 1-D 
problem of §4.3.1. The numerical pulse diffuses too slowly, so that it is slightly too large near the center and 
slightly too small further out. 

Running with the setup described and using the results displayed in Figure 15, we calculate the residual 
as a quantitative measure of the solution quality. With the Euclidean norm, the test exhibits a relative residual 
at the final timestep of 
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0.008. (77) 



4.3.3. Flux Divergence Test 



The flux divergence test was posed by Turner & Stone (2001) as a means of testing the divergence term 
of the diffusion equation. Our problem setup is modeled after that test described in §5.2 of Turner & Stone 
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(2001). An analytic solution is known for the case of the diffusion equation on a unit square in Cartesian 
coordinates with periodic boundary conditions and a constant diffusion coefficient. For this problem, the 
diffusion equation can be solved via separation of variables. Lamb (1995) outlines an analytic solution: 

E e (x,y,t) =2 + e^ 2t {sm(2nx)}{sm(2ny)}. (78) 

For the numerical test, we discretize a unit square in curvilinear coordinates into 100 uniform zones 
in each spatial dimension in a fashion similar to the previous 2-D Gaussian verification problem. The inner 
edge of the grid is located at r$ = 10 4 cm. For the radiation, we use a single energy group, bounded by and 
1 MeV, with the arithmetic average of the group boundaries, 0.5 MeV, assigned as the group-center value. 
Since this is a pure transport test and the radiation microphysical cross-sections are set such that opacities 
are independent of material properties (D = 1), the values of the material variables are irrelevant (other than 
the specification that the medium is hydrostatic). The system is initialized to the differenced form of the 




Fig. 16. — The results of the radiative flux-divergence test showing the maximum relative error in each of 
the spatial zones. In our test calculation, on the final timestep, the maximum relative error that any point 
reaches over the entire domain is about — 1.6 x 10~ 8 . The maximum relative error that any point reaches 
over the entire domain at any time in the course of the simulation is about 1.7 x 10~ 4 . 



analytic solution at a small time, t = 10 s, 

[£ e ]/+(i/2),y + (i/2) =2 + e- 87r2, °sin{27T([x] ; . +(1/2) - M,. +(1/2) )}sin{27T(ty] . +(1/2) - M ; - +(1/2) )} , (79) 

and is left to evolve for 1 second (approximately 100 timesteps). Timesteps are chosen such that At = 
lOOAr/c, where Ar is the spacing of the radial zones. 
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The comparison between the numerical and analytical results is very good. Figure 16 shows the maxi- 
mum relative error of the numerical solution across the spatial domain. This is the maximum relative amount 
by which each spatial point deviates, at any time, during the course of the calculation. Therefore, at a typical 
timestep, the solution is actually better than that plotted. In our benchmark calculation, on the final timestep, 
the maximum relative error that any point reaches over the entire domain is about —1.6 x 10~ 8 . The maxi- 
mum relative error that any point reaches over the entire domain at any time in the course of the simulation 
is about 1.7 x 1(T 4 . 



At the final timestep, using the Euclidean norm, we find our results exhibit a very small relative residual, 

9.1 x 1(T 9 . (80) 
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This number is small at all timesteps throughout the calculation. It sees its maximum value (~ 7 x 10~ 5 ) 
right at the beginning of the calculation. It then declines almost monotonically to the final value shown in 
equation (80). 



4.3.4. Light-Front Propagation 

The light-front propagation problem tests the ability of the algorithm to propagate a front in the opti- 
cally thin (free-streaming) limit where the behavior of the Boltzmann equation is hyperbolic. This presents 
a challenging test for a flux-limited diffusion algorithm. This test was used for the flux-limited diffusion 
algorithm of Turner & Stone (2001) (see their §5.5). 

The problem consists of modeling the propagation of a light front, following it from its start, at the 
inner edge of the domain, until it reaches the halfway point of the domain. The time for this to occur is just 
the straightforward expectation, based on a front traveling at speed c. 

The numerical version of this problem is set up in spherical polar coordinates, at a large radius r , which 
we choose as 10 4 cm. We divide a 1-cm spatial domain in the radial direction into 100 zones, equally spaced. 
A single radiative energy group is also assigned. It is arbitrarily bounded by and 1 MeV and given the 
arithmetic average value, 0.5 MeV, at the group center. We set the total radiative opacity, Kj = 10~ 5 cm -1 , 
which is constant in space and time. Like the other transport verification problems in this subsection, there 
is no emission or absorption. The initial value for the radiation energy density within the computational 
domain is set to a suitably small value (< 10~ 5 erg cm~ 3 MeV -1 ) to avoid encountering possible problems 
with zero or negative energy densities. 

Energy is supplied to the domain by a Dirichlet boundary condition at the inner edge of the domain. 
This is accomplished by fixing all left-hand-boundary zones at all times to 1 erg cm~ 3 MeV -1 . This po- 
sitions the radiation front at x = when t = 0. Since the problem is of limited duration, the boundary 
conditions chosen for x = 1 cm, the right-hand boundary, are largely irrelevant. We run this problem with 
flat boundary conditions (see Appendix H). 
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We also run this test as a nominal two-dimensional problem, with 2 zones in the G direction that span 
10~ 4 rad. Periodic boundary conditions are applied to quantities in this second dimension. 

These initial conditions are slightly different that those used by Turner & Stone (2001), who fix K T 
to achieve an optical depth of 0.01 for the domain. This places our version of the problem in an optically 
thinner medium, since our smaller value of K T gives our domain an optical depth of 10~ 5 . In addition, 
Turner & Stone position the radiation front initially at x = 0.1 cm. 

Our results are shown in Figure 17. This plot corresponds to a time of about 1.7 x 10~ n s. The curve 
shown in Figure \l{a) shows results from our standard version of the test. The timestep is selected as 0.5 
times the minimum transport CFL time, Ar/c, in the problem. Since the zones are equally spaced, this 
corresponds to a timestep of 1.7 x 10~ 13 s. The problem runs 100 timesteps to reach the outcome shown. 
Clearly, the numerical solution is somewhat diffusive. 

Improvement can be obtained when timesteps are reduced. However, even then, the advantage realized 
is limited. This is illustrated in Figures 17(b) and (c), where calculations were made using CFL fractions 
of 0.05 and 0.005, respectively. Although both the radiation fronts are much better resolved than for Fig- 
ure 11(a), 10 and 100 times the amount of computational effort was expended in producing results (b) and 
(c) relative to (a). Furthermore, the resolution improvement in (c) relative to (b) is probably not significant 
enough to warrant 10 times the computational effort. 

In Figure 11(d), we set the CFL fraction back to 0.5, but divide the domain into 1000 zones in the radial 
direction. The combination of more zones, which are thinner and thus decrease the timestep size, results in 
a calculation that requires roughly the same computational work as performed for (c). The resolution of the 
radiation front is also comparable, though the increased spatial resolution noticeably sharpens the front and 
reduces the presence of its numerical precursor. 

An analogous set of tests are performed for light fronts propagating in the X2 direction. These tests 
achieves similar results to those of the x\ direction tests. 

All the results we obtain compare favorably with those of Turner & Stone (2001), who also use the 
Levermore-Pomraning flux limiter in a diffusion scheme to model their radiation. However, in both cases, 
the diffusivity evident in the solutions demonstrates some ultimate limitations in using flux-limited diffusion 
in a regime that is maximally far away from the diffusion limit. 

4.4. Radiation Hydrodynamics Tests 

We have conducted a set of verification tests designed to elucidate the performance of our algorithm 
on radiation-hydrodynamic phenomena. These tests focus on the exchange of energy between matter and 
radiation under a variety of conditions. The first of these verification problems tests the exchange of energy 
in cases where the matter is static. The second tests the exchange of energy under conditions where the 
dynamics of the matter drives the exchange. 
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Fig. 17. — Results of the light-front test in spherical polar coordinates. The coordinate value x = is located 
at radius, r Q = 10 4 cm. The plots show the calculated solutions using various fractions of the minimum 
transport CFL time: (a) 0.5, (b) 0.05, (c) 0.005, and in the case of (d), 0.5, but calculated using 1000 radial 
zones. 



4.4.1. Relaxation to Thermal Equilibrium Tests 



This test, suggested by Turner & Stone (2001), is designed to test the ability of a code to achieve 
thermal equilibrium between matter and radiation by means of emission and absorption. This problem is a 
gray problem where the emission is blackbody in nature. Our multigroup algorithm can accommodate this 
problem by simply setting the number of groups to one and setting the group width to Ae = 1 MeV. The 
gray opacities and emissivity can then be used for this single group to model a gray problem. We consider 
two cases: one where the initial matter temperature is above the equilibrium value and the other where the 
initial matter temperature is below the equilibrium value. In both cases, the initial radiation energy density 
is set to E = 10 12 erg cm~ 3 . The equation of state is described by a y = 5/3 ideal gas with a mean molecular 
weight of /J. = 0.6. The density of the material is taken to be p = 10~ 7 g cm~ 3 and the opacity is taken 
to be K a = 4 x 10~ 8 cm -1 . In the first case, the initial value of the matter energy density is chosen as 
e = 10 2 erg cm~ 3 ; in the second case the initial value is e = 10 9 erg cm~ 3 . The timestep was chosen to be 
At = 10~ n s with the initial time set to to = 10~ 16 s. The numerical results of these tests are displayed in 
Figure 18. 




Fig. 18. — Two trajectories towards equilibrium: (a) via radiative heating and (b) via radiative cooling. In 
both plots, the analytic solution is plotted in red and representative samplings from the numerical solutions 
are indicated by crosses. 



The analytic solution to this problem can be easily obtained by integrating the ordinary differential 
equation defined by the emission-absorption terms: 

de 



where 



cK a E-4nB(T), (81) 

B(T) = j^T 4 , (82) 
um b (v- \)e 
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and cir is the radiation constant and mi, is the baryon mass. If we make the assumption that the radiation 
energy density is approximately constant, then by substituting equation (83) into equation (82), we can 
rewrite (81) as 

f = r-r]V, (84) 

where 77 and j3 are constants. Clearly, the initial value of the matter energy density e , at time to, reaches a 
final equilibrium point <?y = /r\ at some time tf. Thus, we can integrate equation (84) to find 

f' , f e de' 

The indefinite form of the integral on the right-hand-side of equation (85) can be expressed in two forms 
(Selby 1972; Gradshteyn & Ryzhik 1994) 




^"^tan^U ifx^/17 



P/V+xJ Vj3 

(86) 



which depend on the relationship of x to y3 / T7 . Thus, we can rewrite equation (85) to yield an analytic 
relationship between e and t: 

t(e)=t + f(e)-f(e ). (87) 
In this form, it is easy to compute the time t(e) at which a value of the energy density e obtains. 

In Figure 18, we compare the numerical and analytic solutions that we obtain for both the heating 
and cooling cases. In both cases, the solution achieves the equilibrium value, as expected. In the radiative 
cooling case, we note that the numerical results lead the analytic solution slightly for the first few steps. 
This is expected, since the cooling timescale is initially somewhat shorter than the numerical timestep — 
resulting in an slightly inaccurate solution. As time evolution proceeds and the matter radiatively cools 
towards equilibrium, the cooling timescale exceeds the time step, and the numerical solution becomes nearly 
identical to the analytic solution for subsequent times. The cooling-case result makes an interesting contrast 
to the results reported in Figure 3 of Turner & Stone (2001), who indicate that their numerical solution for 
cooling lags the analytic result at early times. The difference between are solution and theirs arises because 
our algorithm does not solve the Lagrangean portion of the gas energy equation simultaneously with the 
radiation-diffusion equation. 



4.4.2. Radiation-damping of Acoustic Hydrodynamic Waves 

This verification test, first carried out by Turner & Stone (2001), is based on the linearized analysis 
of the radiation-hydrodynamic equations carried out by Mihalas & Mihalas (1983) (hereafter MM83). This 
analysis predicts the radiation-damping of acoustically driven hydrodynamic disturbances. The dynamics of 
the acoustically driven wave depend on the ratio of the optical depth of the problem to the damping length. 
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Our test utilizes the same parameters as those employed by Turner and Stone. The gas is an ideal gas 
with adiabatic index of 7 = 5/3. The unperturbed matter density is specified as po = 3.216 x 10~ 9 g cm~ 3 . 
The unperturbed matter internal energy density is given by Eq = 26 020 erg cm~ 3 which corresponds to an 
unperturbed temperature of T = 5.59043 x 10~ 6 MeV and an adiabatic sound speed of c s = 2 998 295 cm s -1 
(10 _4 c). The unperturbed radiation energy density is specified as E e o = 17 340 erg cm~ 3 MeV -1 . We utilize 
a single energy group centered on an energy of 1 MeV with a width of 1 MeV. This choice of grouping 
makes the MGFLD equations behave as though they were gray flux-limited diffusion equations. This choice 
of problem parameters yields a Boltzmann number for the system of 

B ^4yt £ £o = io _ 3 ^ (gg) 
cE e0 

where c s is the adiabatic sound speed of the gas given by 




(89) 

In the linearized analysis of MM83, the system is also characterized by the quantity r given by 

r= — B (90) 



c 



which, for the parameters described above, yields the value of r = 0.1. Finally, the matter is assumed to 
have an absorption opacity given by fc a = 0.4p cm -1 . 

MM83 posit plane-wave solutions to the linearized radiation-hydrodynamic equations which allows us 
to write the density, matter internal energy density, and velocity in the form 

Pl = PoOCpe^'-^, (91) 

Ei=E a E ^ w - kx \ (92) 

0l = v e^ m - kx \ (93) 

where a p , GCe, and Vq are constants specifying the amplitude of the acoustic perturbation. 

Based on these perturbations MM83 derive a dispersion relationship for the time-dependent linear 
radiation-hydrodynamics equations described by the equation 

Az 4 + Bz 2 +C = (94) 

where A, B, and C are complex coefficients given by 
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C=-3zl 



and where we define 



*c s = — 
CO 

and 



(98) 
(99) 



CK a 

?c = • 

CO 

The quantity T Ct is related to the optical depth per acoustic perturbation wavelength. MM83 suggest that 
equation (94) must be solved numerically, but an simple analytic solution of this particular form of quartic 
equation exists. By the substitution u = z 2 equation (94) is reduced to a complex quadratic equation, for 
which the standard quadratic formula applies. The roots are then given by 

z = z r + iZi = ±Vu- (100) 

The wave number k can then be written as 

CO CO CO 

k = -z = -z r - -Zii, (101) 

C<s C*s C*s 

and, thus, the perturbation takes the form 

gi(mt—kx) x/Lgi(mt— 2nxj K) (102) 

where the damping length L and the perturbation wavelength A are defined by 

L=— (103) 
COZi 

and 

A=^. (104) 

COZr 

The complex relationship between the damping length and the perturbation wavelength is illustrated in 
Figure 19 for the aforementioned values of B and r. 

This figure reveals that acoustic waves with a values of T Ct ~ 10~ 3 or x Cs ~ 10 3 should be only slightly 
undamped while those with t Cj ~ 1 should be heavily attenuated. We point out that, despite appearances, 
the two curves do not intersect near x Cs = 10 3 4 . Our curves, however, show a much closer approach at these 
points than do the artist-rendered curves of MM83. 

To formulate test problems based on this analytic solution of the radiation-hydrodynamic equations, it 
is more convenient to work with the variable 

t a = 2%x Cs = fc a A, (105) 



which is the optical depth per wavelength. 
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Fig. 19. — The damping length to perturbation length ratio versus optical depth per perturbation wavelength. 
In (a), the upper (red) curve shows the root corresponding to matter acoustic waves, while the lower (blue) 
solution corresponds to radiation dominated waves. These two curves do not cross near log 10 T Cj ~ 3.4. 
However, they actually approach much more closely than the artist-rendered curves in Figure 4 of MM83 
leads one to believe. An expanded view of the close approach of these curves is shown in (b). 
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We choose three test cases corresponding to values of Ta of 10~ 3 , 1, and 10 3 . The corresponding values 
of co and L are also given in Table 2. By combining equation (102) with equations (91)-(93), and retaining 
only the real portion, we obtain the solutions 

p =p jl + a p e~* /L cos^-^^ j, (106) 

E = E 1 1 + a E e- x/L cos (at - j , (107) 

and 

cot-—). (108) 

For our test problems we choose the amplitude for the perturbations to be one percent in all cases 
giving a p = (Xe = 0Cer = 0.01 and Do = 0.1c s . The problem is spatially discretized by choosing the zoning 
such that the domain is divided into 150 zones, with the zone width of 10 zones per perturbation wavelength 
(Ar = A/10). The left-end of the domain is located at a radius large enough to put the spherical coordinate 
system into the plane -parallel limit. We have retained the same equation of state and energy group structure 
as used in the radiative heating/cooling problem. The timestep is chosen to be approximately a factor of 10 
greater than the heating/cooling timescale in each problem. 

The analytic solution describes a wave moving from left to right across the domain. The density, 
velocity, matter internal energy density, and radiation energy density are initialized to their analytic values 
at t = 0. The perturbations are driven by reseting these quantities to their analytic values in the leftmost ten 
zones, after the completion of every timestep. 

Our numerical results for late-times are displayed in Figures 20-22. The behavior in all cases agrees 
well with the analytic solution. In the optically thick and thin cases, the numerical solution is damped just 
as predicted while the optically translucent case produces a solution that is almost unattenuated — again, as 
predicted. Since the timesteps are much longer than the radiative diffusion timescales, the radiation energy 
density is essentially flat in all cases. If we choose timesteps shorter than the heating/cooling timescales, we 
can also reproduce oscillatory behavior in the radiation-energy density. 

Table 2. PARAMETERS FOR ACOUSTIC WAVE TESTS 





CO 


L 


A 


L/A 


Behavior 3 




(rads- 1 ) 


(cm) 


(cm) 






10 3 


2.423 x 10~ 5 


4.912 x 10 12 


7.774 x 10 11 


6.316 


damped 


1 


2.423 x 10~ 2 


9.775 x 10 11 


7.774 x 10 s 


1257 


undamped 


10~ 3 


2.423 x 10 1 


1.374 x 10 6 


7.774 x 10 5 


1.766 


damped 



The behavior of the material acoustic wave is listed for each set of parameters. 
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Fig. 20. — Acoustic damping for the Ta = 10 3 case. The y-axes represent 100c)X/Xo where 8X = X — Xq 
where X is p, e, v, or E. The region where the solution is reset after every timestep is to the left of the 
vertical dashed line. 
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Fig. 21. — Same as figure 20 but for the Ta = 1 case. 
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Fig. 22. — Same as figure 20 but for the Ta = 10 3 case. 
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4.5. Parallel-Integrity Tests 

The parallel-integrity tests are a set of miscellaneous tests that are intended to confirm that the parallel 
implementation of V2D is correct. As such, they are more in the nature of computer-code implementation 
tests instead of physics tests. Nevertheless, since a portion of this paper is devoted to the parallelization of 
this algorithm, we include a brief description of these tests. The regular execution of these tests is critical to 
ensure that the parallel implementation is bug-free. 

We describe three tests in this section. The first two check for correct parallelization of I/O in the code. 
The third tests for the correctness of parallelization of the entire V2D code on a test neutrino-radiation- 
hydrodynamic problem. There are numerous other parallel integrity tests we run that are not described 
here for reasons of brevity. Although they are described in other sections, we note that the shock tube and 
unit-density sphere tests described in §§4. 1 .2 and 4.2. 1 are performed in such a way that makes them parallel- 
integrity tests as well. These tests are run on different process topology configurations and the results are 
compared to ensure that answers are bitwise identical in all variables regardless of the process topology. 
In fact, many of our verification tests are likewise run on a variety of process topology configurations and 
bitwise comparisons are performed among the results. 

4.5.1. I/O Process Topology Comparison Test 

This test serves as a test of V2D's implementation of parallel I/O, which uses files formatted in HDF5. 
An initial HDF5-formatted file is created via a serial program. V2D is rerun multiple times on diverse 
numbers of processors. Each time, an initial data file is read in via parallel I/O and written back out in 
parallel. Following this, the output files and input files are then compared to verify that they are identical. 

These V2D parallel I/O runs are carried out with several different processor topologies: lxl, 2x2, 
4x4, 8x8, and 16x16, which correspond to 1-, 2-, 4-, 16-, 64-, and 256-processor runs. The processor 
topology refers to the manner in which the contents of each array have been assigned to a processor. In 
each case, the data has been decomposed according the spatial dimensions of the arrays. As an example, 
the 16 x 16-processor run sees a 32x32 array divided into 256 portions — each of dimension 2x2, with one 
portion resident on each processor. The six distinct output files produced by the multiple runs are compared 
to each other. We find that these files are bitwise identical, thus providing confidence in the correctness of 
our parallel implementation of I/O in V2D. 

4.5.2. Parallel Memory-File Comparison Test 

This test is a variation on the test in the previous section. Whereas the purpose of the "I/O Process 
Topology Comparison Test" is the integrity of data when different processor topologies are used to do 
reading and writing, the focus of this test is to confirm integrity of data resident in a file relative to data 
resident in the distributed memory of an executing parallel program. The two tests also stress two different 
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functional processes by which data is used. The first test stresses a characteristic that is important in the post- 
processing of production data, which occurs across diverse platforms and processor counts. This second 
test stresses a characteristic that is important to the checkpointing, suspension, restarting and reading of 
checkpointed data during the process of performing a single production run. 

A parallel program instantiates three four-dimensional arrays of dimension 40x2x 132x 132, which are 
initialized to specific floating-point values and written to an HDF5-formatted file. The file is subsequently 
read back by the same program, and the data read from the file is compared to data resident in memory. 
This test is performed with the same set of processor topologies that are used in the previous test. The 
decomposition of the data is also handled in the same way. In all cases, we find that the data in memory and 
data written to files agree bitwise. 

4.5.3. Neutrino Radiation-Hydrodynamic Process Topology Comparison Test 

The final parallel-integrity test we perform is a test with the complete neutrino radiation-hydrodynamics 
algorithm. This test stresses all parallel aspects of the code — message passing and global operations in the 
hydrodynamic and transport sections of the code, and parallel I/O of a realistic data set. 

The test is simply constructed. For initial conditions, we take a checkpoint file from a supernova 
model. The test is conducted by evolving the model forward in time by a single timestep on a number 
of different processor topologies using a deterministic global reduction operation for global summations. 
As we discussed in the gravitational unit-sphere-verification-test section, the standard MPI_ALLREDUCE 
global reduction summation operation is non-deterministic and can yield slightly different answers with 
each run. Therefore, for testing purposes, we use a deterministic global summation operation. We wish to 
emphasize that the use of MPI_ALLREDUCE for calculating a global minimum CFL time is deterministic 
and there are no run-to-run variations in this quantity. For the test, a final output file is produced by each run 
with varying process topologies. The test is deemed successful if all variables in the checkpoint files are in 
exact bitwise agreement. 

The deterministic global summation operation that we employ uses a brute force approach in which 
all of the values involved in the summation are systematically communicated to process where they are 
summed in the same order that they would be summed if the code were purely sequential. Although this 
approach is slow, it guarantees that the sum is identical to all bits no matter what process topology is 
employed. 

Although impractically slow for production work, deterministic global summations are needed for in- 
tegrity tests to eliminate the possibility of message passing errors. With the deterministic global sum in 
place, the output files produced by each processor topology of this test should be bitwise identical. Any de- 
viation indicates a programming error in the parallel implementation of the algorithm. Hence, deterministic 
programming is an invaluable development tool in achieving and maintaining consistency and correctness 
of code. 
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In performing this test, we set up the problem in spherical coordinates with 256 radial zones, 32 angular 
zones, 20 neutrino groups for each flavor of neutrino and corresponding antineutrino. V2D is then run for 
a single timestep using a number of processor topologies: lxl, 2x2, 4x4, 16x4, 32x8, where the first 
number indicates the number of processors across which the radial portion of the domain is divided, and the 
second number is the same for the angular portion. A standard V2D checkpoint file is written at the end of 
each of these runs and compared. With the deterministic global reduction in use, we find that all result files 
are bitwise identical. 



4.6. Solver Tests 

The Newton-Krylov solver algorithms used to solve the implicitly discretized radiation-diffusion equa- 
tions are also subjected to regular verification and unit testing. We construct nonlinear and linear test prob- 
lems that are mathematically simple and for which solutions are known. The solvers are then tested on these 
problems and the numerical solution is compared to the known solution to ensure that they agree within ac- 
ceptable tolerances. These tests are carried out on a single processor runs and with multiple processors. The 
solutions are compared with the known solutions in both cases in order to ensure that no parallel implemen- 
tation errors are present. We also have carried out numerous tests using both MPI_ALLREDUCE calls and 
our own deterministic global summation to ensure that the slight process-to-process variability present in 
the summation of the MPI_ALLREDUCE poses no problem for the solution of the implicit systems arising 
from the implicitly discretized radiation-diffusion equations. We have never found any physically significant 
difference in result between the two methods for global summations. 

5. Conclusions 

We have developed a numerical algorithm for 2-D multigroup neutrino-radiation-hydrodynamics that 
is especially suited to modeling stellar core collapse, core collapse supernovae, and proto-neutron star 
convection. However, this algorithm is sufficiently general as to be useful for a variety of 2-D radiation- 
hydrodynamic applications involving radiation other than neutrinos. Our algorithm solves the combined 
set of comoving-frame neutrino-radiation hydrodynamics equations, where radiation is treated in a multi- 
group flux-limited diffusion approximation. The combined set of equations is evolved in time using an 
operator-splitting approach. The nonlinear diffusion equations are evolved implicitly in pair-coupled fash- 
ion while hydrodynamic portions of the operator-split equations are evolved explicitly. The explicit portion 
of this approach draws on a numerical scheme developed by Stone & Norman (1992a) and Stone, Mihalas, 
& Norman (1992). The implicit differencing of this scheme is motivated by the work of Turner & Stone 
(2001). However, the multigroup treatment of the radiation distribution, the pair-coupling between particles 
and antiparticles, the Pauli blocking effects, the solution of implicitly differenced diffusion equations via 
Newton-Krylov iteration, and the means of including an arbitrary equation of state have not been heretofore 
addressed in such an algorithm. 
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The numerical scheme has been described in sufficient detail to allow a complete implementation of the 
algorithm by the reader. We have discussed it's parallel implementation on distributed memory architectures 
via spatial domain decomposition. This includes message-passing strategies to enhance the scalability of the 
algorithm. 

In order to benchmark the numerical performance of the algorithm, we have presented the results 
of a suite of important verification tests. These tests stress diverse portions of the algorithm and include 
hydrodynamic, radiation-transport, radiation-hydrodynamic, gravitational, and parallel integrity tests. Our 
objective has been to establish the numerical accuracy of our algorithm both as a whole and by individual 
aspects. This verification process is vital to establish the integrity of future simulations. This process also 
helps insure that simulations carried out with this algorithm can be properly evaluated by other researchers. 
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A. Finite Differencing Notation 

All discretized quantities are denoted by the use of square brackets around the symbol for the quantity. 
The location in time at which quantities are defined are denoted through the use of standard superscript 
notation and is described in §3.4. The location at which a discretized quantity is located in space or energy 
is denoted through the use of standard subscript notation as follows. Any subscripts inside the square 
bracket denote a component of a vector, tensor, metric coefficient, or a spectral (radiation-energy-dependent) 
quantity. energies.Any leading superscripts inside a square bracket are used to label the specific flavor of 
neutrino-related quantities, e.g., electron, muon, or tauon flavors. 

A coordinate, or a quantity derived from a coordinate, which we generalize by £ , is naturally defined 
either at cell centers or at cell edges in each dimension. We enumerate cell edges with integer values and 
cell centers with half-integer values in each dimension. With this convention, cell center i + \ is surrounded 
by cell edges i and i + 1 . Figure 1 illustrates the 2-D spatial mesh, including the enumeration of cell edges 
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and cell centers in each dimension. The spectral or radiation dimension uses an identical enumeration (see 
Figure 2). 

Discretized coordinates are denoted by a single subscript. Therefore, our notation for the discrete 
analog of coordinate E, at a cell edge is 

E, at the ith cell edge ->[£],■. (Al) 

The value of E, at a cell center is given by 

£ at cell center i + 1 [£],- +(1/2) . (A2) 

In these cases, E, might be spatial coordinates, which we generically refer to as x\ and x 2 ; an energy-space 
coordinate, e; differences of coordinates, metric coefficients, g2, gn, and #32 or their spatial derivatives; 
or the area and volume elements, AA and AV, respectively. We use the terms "group" or "energy group" 
synonymously to refer to cells in the spectral dimension. Usually when we use the terms "cell" or "zone," 
we will be referring to the spatial mesh. 

In general, scalar quantities are defined at cell centers. A generic scalar quantity y(x\,X2), that is a 
function of spatial coordinates x\ and x%, has a discretized analog that is notated as 

Y at the cell center (i+ \,j + \) -► Mi+(i/2)j+(i/2) • ( A3 ) 

Figure 1 illustrates the locations for such cell-centered discretized variables. If the scalar quantity is spectral 
{i.e., a function of the spectral energy e), such as the spectral radiation energy density, we denote it by y e . 
The discretized analog of this quantity, when cell centered in the spatial and energy dimensions, is denoted 

by 

Ye at energy group k+ I and cell-center (i + \j+ \) [Ve] t +(i/2),i+(i/2)j+(i/2) ■ ( A4 ) 
Our convention is that spatial indices follow any spectral indices on all discretized quantities. Figure 2 
shows the location at which such quantities are defined in the 3-D (two spatial dimensions and one spectral 
dimension) mesh. 

A general vector quantity, a, is decomposed into its two components, 0\ and a%. These quantities are 
usually defined on the respective cell faces in the direction of the component. Thus, we adopt the following 
notation: 

<7i at cell face (ij + \)^ [oi] w+ (i/ 2 ) . (A5) 
o 2 at cell face -> [o2] i+(1/2 )j • (A6) 

Note that indices inside the square brackets refer to components of the vector. Figure 1 depicts the locations 
at which the discrete analogs of the vector components are defined. 

If the vector quantity is a function of spectral energy {e.g., the spectral radiation flux density), which 
we denote in general by <7 e , and its components by a e j and a e ,2 , then it is face-centered in the spatial 
dimensions, but group-centered {i.e., cell-centered) in the energy dimension. Thus, 

<T e ,i at group center k+ X - and cell face -> [ cr e,i]fe+(i/2),ij+(i/2) » ( A7 ) 
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a £; 2 at group center k + i and cell face -> [ (J e,2] (t+ (i/ 2 ),i + (i/2)j • ( A8 ) 

Figure 2 illustrates the location at which such quantities are defined. 

Most derivatives are calculated via finite differences, and these places will be made obvious throughout 
the text that follows. For some quantities, however, we calculate derivatives analytically and use these ana- 
lytic expressions in the code. In the present paper, these analytic expressions are often used for derivatives 
of the metric scale factors (see §3.2). In differenced expressions, such derivatives are notated as shown in 
the following example 

when evaluated at grid point i — ► 

ax\ 



(A9) 



B. Discretized Variables 

In this appendix, we present a partial list of the more important physical quantities and the discretized 
notation for these quantities. The conventions for the finite-difference notation we utilize is discussed in 
detail in Appendix A. The physical quantities and their discretized analogs are listed in Table 3. Figures 1 
and 2 show the locations at which these variables are defined. 



C. Generalized Coordinates, Areas, and Volumes 

In this section, we describe how the covariant formulation introduced in §3.2 is implemented in our 
algorithm. In adopting the formulation of Stone & Norman (1992a) for generalized coordinates, we also 
adopt a portion of their notation. Although there are other possible coordinate systems that could be used 
with this method, we focus here on the three most common systems: (i) Cartesian, (ii) cylindrical, and (iii) 
spherical. 

Paralleling Mihalas & Mihalas (1984) and Stone & Norman (1992a), we write the metric in a general- 
ized orthogonal coordinate system as 

ds 2 = { gl ) 2 dxj + (g 2 ) 2 dxl + (gsignfdxj. (CI) 

For Cartesian, cylindrical, and spherical coordinates, Table 4 gives the values of these quantities and their 
derivatives. In discretized form, these quantities are used both on both the integer and staggered, half-integer 
mesh. 

For Cartesian coordinates (x,y), all metric coefficients are unity: 

[gl]i=[g2\i = [g3l]i = [g32]j=l-, 
Mi+(l/2) = [#2];+(l/2) = [g3l]/+(i/ 2 ) = [g32]y+(i/ 2 ) = 1 (C2) 

and all derivatives are obviously zero. 
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Table 3. Sample Discretized Quantities 



Physical Quantity 


Algebraic 
Symbol 


Quantity Type 


Notation 


spatial coordinate: cell edge 


X\ 


independent variable 


[*iL 


spatial coordinate: cell center 


X\ 


independent variable 


1 i \ / / 


spatial coordinate: cell edge 


X2 


independent variable 


[X2]i 
'J 


spatial coordinate: cell center 


X2 


independent variable 


[X2l ;+Cl I7\ 


spectral energy coordinate: 


e 


independent variable 




group edge 








spectral energy coordinate: 


e 


independent variable 


l £ ]k+{l/2) 


group center 








g2 metric coefficient: 


82 


function of independent 


[82], 


cell edge 




variable, x\ 




g2 metric coefficient: 


82 


function of independent 


\S2\ i+ (l/2) 


cell center 




variable, x\ 


mass density 


P 


cell-centered, intensive 


[p]i+(l/2)J+(l/2) 


electron fraction 


Y e 


cell-centered, intensive 


[ y e]i+(l/2)j+(l/2) 


electron number density 


n e 


cell-centered, intensive 


M;+(l/2)j+(l/2) 


matter internal energy density 


E 


cell-centered, intensive 


[ £ ]i+(l/2)J+(l/2) 


matter pressure 


P 


cell-centered, intensive 


[ P ]i+(l/2),;+(l/2) 


matter temperature 


T 


cell-centered, intensive 


[ r ]i+(l/2)j+(l/2) 


X\ comp.: matter velocity 




face-centered, vector 


[ U l]«J+(l/2) 


X2 comp.: matter velocity 


v 2 


face-centered, vector 


M;+(i/2)j 


x\ comp.: matter momentum 


Si 


face-centered, vector 


[ s l]iJ+(l/2) 


density 








X2 comp.: matter momentum 


S2 


face-centered, vector 


h]i + (l/ 2 )J 


density 








spectral radiation energy density 


E e 


cell-centered, intensive 


[ £ £]*:+(l/2),i+(l/2)j+(l/2) 


x\ comp.: spectral radiation 


Fe,l 


face-centered, vector 


[ F e-4+(l/2),i,;+(l/2) 


flux density 








X2 comp.: spectral radiation 


Fe,2 


face-centered, vector 


i F £2] k+( i/ 2 ),i+(l/2)J 


flux density 








1,2 entry: spectral radiation 


{PJl2 


cell-centered, tensor 


[{ P e}l2] fc+ (i/2), i+ (l/2)J+(l/2) 



pressure 



For cylindrical coordinates (z,9t), the metric coefficients are 



\gl\i=\g2]i = \g3l\i = h 
kl](+(l/2) = [g2\i+(l/2) = [S3l]/+(i/2) = !, 

fe 2 ]y = [91]; ; [g32]y + (i/ 2 ) = P*];+(l/2) ■ 



with non-zero derivatives 



dg32 




dg32 


dx\ 




dxi 



1. 



Ji+(l/2) 



(C3) 
(C4) 



Similarly, for spherical coordinates (r,9), on the integer mesh, we have 

M; = l; [82] i = [r] i ; b] ; = M ; ; [g 3 2] j = sm[e] j , (C5) 
while on the half-integer mesh, we have 

Hi+(l/2) = ^ b]i + (l/ 2 ) = Mi+(l/2) ; bl]; + (i/ 2 ) = Mi+(l/2) ; [gv]j+(l/2) = sin[0] y+(1/2) . (C6) 

The non-zero derivatives, when differenced, are 

r) Pi 1 1 r a Pi 1 

= 1 (C7) 



~dg2~ 












dg3l 


dxi _ 




dx\ 


'•+(1/2) 


dx\ 


i 


dx\ 



-I '+(1/2) 



and 



32 



dx? 



:cos[0] ; .; 



J J 



dg 



32 



dxi 



7+(l/2) 



COS [0] 



;'+(i/2) • 



(C8) 



We also make frequent use of the 2-D spatial volume element, dV, which in general can be decomposed 
such that 

dV = dT 1 dT 2 T 3 . (C9) 
The decomposed volume elements, the JT's, are given by 

dYi = g 2 g3i dx\ ; <iT 2 = g32 dx 2 . (CIO) 

The value of T3 depends on the specific coordinate system. 

In Cartesian coordinates, the decomposed volume elements and their discretization are trivial: 

dT x = dx; dT y = dy; T z = 1. (Cll) 

When differenced, they are defined in the code as 

[ATj,- = [x] i+(1/2) -[x] ; ._ (1/2) (C12) 
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Table 4. Coordinates and Metric Scale Factors for Selected Coordinate Systems 



Quantity 


Symbol 


Cartesian 


Cylindrical 


Spherical 


coordinates 


X\ 


X 


z 


r 




x 2 


y 




e 






z 


# 




metric scale factors 


Si 


l 


1 


i 




82 


l 


1 


r 




#31 


l 


1 


r 




#32 


l 


91 


sin(9 


scale factor derivatives 










1 




dg-n/dxi 








1 




dgii/dxi 













dg 2 /dx 2 













dg3i/dx 2 













dg3l/dx 2 





1 


cos B 


decomposed volume elements 


dTi 


dx 


dz 


r 2 dr 




dT 2 


dy 


9td9t 


sinddd 




T 3 


1 


2% 


2n 


area elements 


dA\ 


dy 




2nr 2 sin d dd 




dA 2 


dx 


2n3id z 


2nrsind dr 



-74- 



[AV] 




i ;+(l/2) i+1 



Fig. 23. — Location of the volume element that corresponds to the volume of the computational cell 
[AV],- + (i/2)j+(i/2)- Volume is defined as a cell-centered quantity and is therefore defined at the junction 
of the half integer mesh at the center of a cell. The third dimension, which makes this enclosure a volume, 
is orthogonal to the page and is in the direction of the assumed symmetry. 



and 

[AU- = [>W)-MMi/2)- (C13) 
Cell-centered and cell-edge coordinates in the x direction are usually related by 

M.-+(i/2) = ^{Wm + Mi}, (C14) 
with the analogous relationship for the y coordinates. Cell-centered volume elements are also used: 

[ATj i+(1/2) = H I . +1 -[x] j (C15) 

and 

[AT y ] i+(1/2) = M m -M,- ( C16 ) 
In cylindrical coordinates, with azimuthal symmetry, these elements become 

dT z = dz; dYx=Kdft; Y#=2n, (C17) 
which, upon discretization, become 

[^z] i =[z] i+{1/2) -[z] i . {1/2) , (C18) 

[ATsdj = 2% (pt] y+(1/2) - Pt],-- ( i/2)) • (C19) 
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Cell-centered and cell-edge coordinates in the z direction are usually related by 

[z]/+(i/2) = ^{[z]m + [z]i}> ( C2 °) 
with the analogous relationship for the 9^ coordinates. 

In spherical coordinates, with azimuthal symmetry imposed, these elements become 

dT r = r 2 dr; dT e =smddd; T^=2n. (C21) 
When differenced, they appear in the code as 

[AT,], = I {(M i+(1/2) ) 2 + [r],. +(1/2) [r] M1/2) + (M f _ (1/2) ) 2 } ([r] i+(1/2) - [r] i _ (1/2) ), (C22) 

[AT e ] . = - cos [0] . +(1/2) + cos [0]j_ {1/2) . (C23) 

We note that writing the volume element in the form used in equation (C22) avoids the computation of the 
difference of two cubes, which can lead to a substantial loss of precision. The choice of relationship between 
cell-centered and cell-edge coordinates in spherical geometry is quite important. Following Monchmeyer & 
Miiller (1989) we usually choose 

M*(i/2) = ^JlW^y + irUMi + drm (C24) 
for radial coordinates, while choosing 

to relate cell-centered and cell-edge coordinates for the angular dimension. Cell-centered volume elements 
are also used: 

[ATr]^ = 3 {(Mm) 2 + Mm Mi + (Mi) 2 } (Mm - M,-)> (C26) 
[ AT e]y+(i/2) = - cos [0] j+l + cos [6] j . (C27) 

Also of importance are the components of the differential area vector, dA, in the x\ and x 2 directions. 
These are given, in general, by 

dA\ = Y 3 g 2 g 3 i dY 2 ; dA 2 = T 3 g 3 ig 32 dxi. (C28) 
Since these elements are components of a vector, these quantities become face-centered when differenced 

[AAi] iJ+(1/2) = T 3 Mi fell, [ A H+0/2) (C29) 

and 

l^iMWJ =T 3 bl]i+(l/2) k32],-([xi] i+1 " [*lL0- (C30) 
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. [AA2];+(i/2)j+i 



7+1 



[AA,] 4 



;+(l/2) 



[AAi] 



i+l,j+(l/2) 



i+l 



[AA 



2ji+(l/2)J 



Fig. 24. — Location and orientation of the components of the differential area vector with respect to the 
volume element they enclose. The third dimension, which makes these cell "faces" areas, is orthogonal 
to the page and is in the direction of an assumed symmetry. Area, being a vector quantity, is defined at 
cell faces. In our staggered mesh algorithm, this means that vector components originate from different 
locations. 
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Figure 24 shows the location and orientation of the components of these area vectors with respect to the cell 
volume they enclose. 

For Cartesian coordinates, the area elements and their discretization are again trivial. 

dA x = dy; dA y = dx, (C31) 

which, when differenced, give 

[AAJy+(i/2) = M;+i " Mj (C32) 

and 

[ AA y]/ + (i/2)j = Hi+i-H,- ( C33 ) 
In cylindrical coordinates, imposing symmetry in the azimuthal direction, the area elements are 

dA z = 2jflHd3t; dA % = 2nftdz, (C34) 

which, when differenced, give 

[AAj,- J+(1/2) = 2n m j+{l/2) - Pfy) (C35) 

and 

[AA*] i+(W = 2tt ([ Z ] i+1 - [ Z ],.). (C36) 
In spherical coordinates, imposing azimuthal symmetry, the area elements are 

dA r = 2nr 2 sin dA e = 2nrsm ddr, (C37) 

which, when differenced, give 

[AM-,y+(i/2) = 27r([r] ! .) 2 (-cos [0] ; . +1 +cos [0] y ) (C38) 

and 

[AA e ] /+(1/2)J =2^[r] I . +(1/2) sin[0] ; .([r] i . +1 - [r],.). (C39) 

D. Discretization of Advection Equations 

In this appendix we discuss the solution of equations (32), (34), (37), (42), and (48). Our numerical 
method for these equations is that of Stone & Norman (1992a) and the first subsection this appendix merely 
restates their approach for the convenience of the reader. In later subsections we write down the update 
equations in a convenient form for numerical implementation. 
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D.l. The ZEUS Advection Scheme 



Since equations (32), (34), (37), (42), and (48). have the similar form, we can generalize these expres- 
sions. First, for advection of a scalar quantity, we write 



~di~ 



+ V-(yv) = 0, 



(Dl) 



advection 

where iff represents the scalar field being advected. Similarly, for advection of a vector field a, we write 



do 

L 01 U advection 



+ V-(<Tv) =0. 



(D2) 



In solving the advection equation, it is common to solve finite-difference approximations to integral forms 
of equations (Dl) and (D2). For advection of a scalar field, this takes the form 



and for advection of a vector field, 



d_ 

dt 



dt 



\ ydV = - yv-dA, 
Jv Ja 



I adV = - [ G\-dA, 

Jv Ja 



(D3) 



(D4) 



where V is an arbitrary volume whose bounding surface is A. 

Following Stone & Norman (1992a) we directionally split the update of all equations of the form of 
(D3) and (D4) into separate sequential updates of the variables due to advection in the x\ and x-i directions. 
We alternate the directional order of advection updates (x\ updates followed by X2 updates and vice versa) 
with each timestep. A complete directional "sweep" is performed by solving all of the advection equations 
corresponding to one direction, followed by a second sweep solving the equations corresponding to the 
orthogonal direction. 

The consistent advection scheme of Stone & Norman (1992a) ties the advection of all variables to the 
mass flux. Therefore, we first consider the update to the continuity equation due to advection in the x\ 
direction. We can write the discretized version of equation (32) as 



[AV] 



i+ (l/2)J+(l/2) 

At 



[p]"+a/2)j+(i/2) [pr 



i+(l/2)j+(l/2) 



" (^(P)]" + + "y + (i/2) [AAi] mj+(1/2) - [*i(p)]J# (1/2) |AAi] ; j+(i/2)) , (D5) 

where the bounding volume is now [AV] i+ (i /2)j+(i/2)> tne volume of an arbitrary cell in the computational 
mesh. The AA,'s are the areas of the cell faces that, collectively, bound AV. (These cell-face areas are 
defined in Appendix C.) The quantity [^i(p)],-j+(i/2) represents the mass flux density (or mass flux per unit 
area) in the x\ direction across the cell face at (/, j + (1 /2)). The superscripts indicate the level of temporal 
update, with n + a being the time-level at the beginning of the x\ advection update, and n + j3 being the 
time-level upon completion of the update. If the x\ update is carried out before the X2 update we have a = 



j+1 



Ml/2) 



,/+(l/2) 



i+1 



7-K1/2) 



i+lJ+(l/2) 



[^W]'+(l«)jtl 



[F 2 (V|/)] i+( i/2,j 



i+(l/2) i+1 



(+(1/2) i+1 



Fig. 25. — Definition of locations for the x\ and X2 components the flux densities used in the advection of 
a general, cell-centered, scalar quantity [v]i+(i/2)j+(i/2)' The shaded region indicates the volume element, 
[AV]i + n/2)j+(i/2)> use d in the differenced integral advection equation (eq. [Dl]). Since flux density is a 
vector quantity, its components are defined at the four face centers bounding the volume element. 
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and P = a; if the x 2 update is carried out first we have a = a and j3 = b. The update equation for mass 
advection in the x 2 direction is given by 



[ A H-+(l/2)J+(l/2) f.^n+S ln] n+ 7 \_ 

A{ ^PJ/+(l/2)J+(l/2) iPJi+(l/2)j+(l/2)J - 

" (^(P)]^ /2)J+ i [AA 2 ] ; - +(1/2)J+1 - [F2(p)t{ l/2)J [AA 2 ],. +(1/2)J ) , 



(D6) 



where [i 7 2(p)]i+(i/2)j represents the mass flux density x 2 direction across the cell face at (/ + (1/2), j). As 
in the xi update the superscripts indicate the time level at the beginning and end of the update: if the x\ 
update is carried out before the x 2 update we have 7 = a and 8 = b; if the x 2 update is carried out first we 
have 7 = and 8 =a. 

Readers who are contrasting our equations with those of Stone & Norman (1992a) will notice that 
the ZEUS algorithm time-centers many of the quantities that appear in the advection equations. This is 
necessitated by the possibility that the grid may be moving. Since we are considering only grids that are 
fixed in time, the advection prescription here is more straightforward. As a consequence of fixing the grid, 
the spatial coordinates, the volume and area elements, as well as the metric coefficients, are all constant in 
time. 

Simple algebra allows equations (D5) and (D6) are readily solved for the new values of the density in 
terms of the old values and the mass fluxes which are computed base on those old values. 



The mass fluxes are constructed as 

[*l(p)]y+(l/2) = [^l(P)k;+(l/2) Ny+(l/2) 

and 

[Hp)]i+{l/2),j = [^2(p)L- +( i /2 )j[H'+(l/2)J> 

where the ^'s are van Leer mono tonic upwind advection functions (van Leer 1977) given by 

[ U lk;+(1/2)M 



[My)] 



i,;+(l/2) 



Mi-(l/2)J+(l/2) 
[V];+(l/2)j+(l/2) 



+ [Si(V)];-(i/2)j+(i/2) 

- [Sl(V)];+(l/2)j+(l/2) 



1 + 



[*lJ;-L*lJi-l 

: J+(l/2)' 



if [ui 



U+(l/2) 



(D7) 
(D8) 

>0 



if [ u iky+(i/2) < 0, 
(D9) 



where 



[Si(v)L- 



+(l/2),;+(l/2) 



[AvL 



;+(i/2) 



i,y+(i/2) 



[VL'+(3/2)J+(l/2) - [^]/-(l/2)J+(l/2) 

if [ A Hj + (i/2) 



[Ayr] 



i+lJ+(l/2) 



>0 



(D10) 







otherwise, 
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and 

[^+(1/2) = [V] ; - + (l/2)J+(l/2) " M H l/2).; + (l/2) • ( D11 ) 

Here, and in equations that follow, we have temporarily suspended temporal superscripts whenever writing 
expressions at a generic time-level. Equations for updates will retain the appropriate time superscripts. For 
reasons of brevity, we will not write down the analogous expressions for the x 2 direction. However, Table 5 
gives a complete prescription for constructing ^(y) from the J^i(y) definition. 

The advection equations for scalar quantities other than the density are discretized as 

[ Ay ]i+(i/2),y+(i 



At 



(1/2) f\vr] n+P -\vf\ n+a \- 

ijn-+(i/2)j+(i/2) m-+(i/2)j+(i/2) ) - 



- ([Mv0]£j+(i/2) [AAi] /+1J+(1/2) - [FimTiHim Ni] u+(1/2) ) , (D12) 

for the x\ update and 



[ A H-+(l/2)J+(l/2) / ]n+5 u ] n+y \ _ 

X t ljn-+(l/2)J+(l/2) " lV\i+(l/2),j+(l/2)) ~ 

- (fc(vO]?(Wi [ AA 2] I - + (i/2) J+ i - mw)] tlmj V^Mmj) (D13) 



for the x 2 update. 

One of the hallmarks of the ZEUS advection strategy is the use of 'Norman's consistent advection 
scheme (Norman et al. 1980), which ties the fluxes for all advected quantities to the mass flux. For scalar 
quantities other than the density, the fluxes are written as 

p. 



[Mv)L-j +( i/2) 

and 



[*i(p)]ij +( i/2) (° 14 ) 

U+(l/2) 



■y 2 



(- 

VP 



i+(l/2)J 

When equations (D14) and (D15) are substituted into equations (D12) and (D13) we can solve for new 
values of all scalar quantities in terms of old values and the fluxes which are calculate in terms of those old 
values. 

Advection of vector quantities proceeds slightly differently because of the location where vector com- 
ponents are defined on the staggered mesh. The analog of equations (D5) and (D6) for the advection of Ci , 
the x\ component of vector a, in both the x\ and x 2 directions is 

- ([Fl(Ol)]'i+"/2),j+(l/2) [ AA l]/+(l/2),y+(l/2) - [ F l( a l)]HV2)J+(l/2) [ AA Ai-{ll2),j+(l/2)) ( D16 ) 



and 

[Ay 

At 



\- AV h+(l/2) f r ]n +S r f+r \ _ 

([^(oi)]^ [AA 2 ],, +1 - [^(oi)]^ [M 2 ]y) } , (D17) 
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Fig. 26. — Locations of the flux densities needed to perform an x\ advection sweep of a vector quantity a. 
The fluxes used to advect component Oi are shown in (a), while the fluxes for component 02 are shown in 
(b). The shaded regions indicates the volume elements used in this x\ sweep. Volume element, [AV],- ; + (i/2), 
shown in (a) is used in the Ci updates, while in (b), [AV] i+ ^/ 2 ),j i s use< i i n tne °2 updates. These are 
different volume elements than used in the advection of a scalar. This is because of the face-centering of 
vector quantities, whose component locations are indicated by blue arrows. Because the volume elements 
are different, the locations of the flux densities are also different. The interpolations used to obtain them 
are given by equations (D18), (D19), and their 02 analogues. The face-centered locations that are used to 
calculate the interpolations are shown as red circles in the figure. 
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Fig. 27. — Locations of the flux densities needed to perform an %2 advection sweep of a vector quantity O. 
The fluxes used to advect component C\ are shown in (a), while the fluxes for component a% are shown 
in (b). The comments made in the caption to Figure 26 regarding volume elements, the locations of flux 
density vectors, and fluxes used in interpolation formulae also apply here. 
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where the fluxes are given by 

[ F l( a l)]i+(l/2)J+(l/2) 

and 



0\ 

P /J/+(l/2)J+(l/2) 
<?1 



[Flip)] 



i+(i/2),y+(i/2) 



[F2(P)] Lj . 



(D18) 



(D19) 



We note that fluxes in the x\ direction (the Fi's) are here needed at cell centers and the fluxes in the X2 
direction (the F 2 's) are needed at zone vertices — locations where neither of these quantities is naturally 
defined. (Figures 26(a) and 27(a) show the location of fluxes needed to the advect a\.) We use the following 
averages to supply the needed values. 



and 



[*i(p)]i+(i/2),/f(i/2) = \ ([^i(P)]y+(i/2) + ^i(P)]/+ij+(i/2)) 
\UP\j = \ ([^(P)] ! - +(1/2)J + [F 2 (p)],_ (1/2)J ) . 



(D20) 



(D21) 



The van Leer interpolation function for advection of the x\ component of a vector in the x\ direction, 
Jx[c\),\s, given by 



lj+(l/2)J+(l/2) 



[^i]«j+(i/2) + [5i(<yi)]yH 



(1/2) 



^ ([m],-J + (l/2) + [m],- + lJ + (l/2))^ ' 
2 (t x l]i+(l/2) - [ x l]j-(l/2)) 

if [wi]y+(i/2) + [wi]i+i,y+(i/2) > 



(^i] I -,y + (i/2) + [ui],- + ij +( i/2))^^ 

[^lJ/+lJ+(l/2) - [Ol(ffl)Ji+lJ+(l/2) I 1 + -7— — \ 

2 ^l*lJi+(3/2) - l*lJi+(l/2) J 



if [ u i]«j+(i/2) + [«i]/+i,y+(i/2) < 0, 



(D22) 



where 



[ A<J 1 ] i- ( 1/2) ,j+ ( 1 /2) t Acr l ] i+ ( 1 /2) J+( 1/2) 



[5i(o-i)] ( . J+(1/2) 



Nmj+aA) - [ffl],--i,y+(i/2) 

if [A(7i] I -_ (1/2)J+(1/2) [A<7i] ( - +(1/2)j . +(1/2) > (D23) 

otherwise, 



and 



[ AcJ l]i+(l/2)J+(l/2) - [ ff l]i+l,y+(l/2) - [ a lk./+(l/2) 



(D24) 
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The interpolation function needed for the calculation of the fluxes in the x 2 direction can be obtained by the 
substitutions described in Table 5. 

The analogs of equations (D16) and (D17) for the x 2 components of a vector are straightforward to 
obtain: 



[ Ay ]i+(1/2)J (r ]n +S 



At 



(N?-KV2),y-N"-Hl/2)j) 



(D25) 



and 



- (^(^2)]^ /2)J+(1/2) [AA 2 ] ( . +(1/2)J+(1/2) - [F 2 (<r 2 )]^ /2)! ._ (1/2) [AA 2 ] ; . +(1/2)J _ (1/2) ) . (D26) 
For details we direct the reader to Stone & Norman (1992a). The fluxes are given by 



02 



[*i(p)]y + [*i(p)]mj 



(D27) 



and 



[/ r 2(a 2 )],- +(1/2)J+(1/2) 



P /J/+(l/2)J+(l/2) 



([^(P)] i+ (i /2 ),y + i + [^(P)] i+ (i/ 2 )j) • (D28) 



The van Leer interpolation function needed for the calculation of the fluxes in the x\ direction, ^\ (a 2 ), 
is given by 



i(^)]f+(i 



i+(l/2)j+(l/2) 



Ni-(l/2)J+[5l(^)]i-(l/2),y 



' t ([m],-,j- ( i/2) + [m],-, 



;'+(i/2) 



A? 



2 ([ X l]«+(l/2)-[^l] I --(l/2)) 

if [wi]y-(i/2) + [wi]y + (i/2)>0 



Ni+(l/2)J-[5l(^)]i+(l/2),; ! + 



i ([<Mi/2) + N,-,j + (i/ 2) )Ar 



2 (t X l]i+(3/2) ~ [ JC l]/+(l/2)) 

if [wi]y-(i/2) + [«i]w + (i/2)<0, 



(D29) 
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where 



[5l(^2)] /J+(1/2) 



= < 



[Ag 2 ],--ij[Ag2],-,j 

[<*]i+(l/2),j-fa]i-(3/2),j 



if [A<y 2 ] I ._ 1J [A(y 2 ] I . y >0 



otherwise, 



(D30) 



and 



[Aa 2 ], 7 . = [o 2 ] iHl/2)J - [<72] H i/ 2)J • (D31) 
The analogous ^2 function can be obtained by the appropriate substitutions as described in Table 5. 



D.2. Update Equations for an x\ Advection Sweep 

Using the discretization scheme from the previous subsection and substituting in the appropriate for- 
mulae for the cell interface areas, we arrive at a set of update equations for a x\ sweep. 



Density update: 



[P]"+(1/2)j+(1/2) ~ [P]"+(l/2)J+(l/2) " At 



Wi+i [g^h+i [^i(P)]"+ij+(i/2) - [Si\i [g3i]i [Fi{p)]"j+ { 



/+(l/2) 



[ATi],. +(1 



Electron flux: 



Fi(*/p)Eri 



U+(l/2) 



P 



Matter-internal-energy-density flux: 



[Fi(£/P)E(1 



7+0/2) 



1 n+a 

<J+(l/2) 

n+a 
-I U+(l/2) 



'+(1/2) 



^(P)]"Ml/2) 



xi -momentum-density flux: 

^(^/P)]- + + (l/2) J+(1 /2) = (^)]]^ / 2) J+(1 /2) ^(P)]]^/2) J+(1 /2) 

x 2 -momentum-density flux: 



^(^/p)]^=[^( u 2)]-r^(p)]"r 



Radiation-energy-density flux: 

[ F d E e/p)] n kt(l/2),i,j+(l/2) 



Si 



-i n+a 



(D32) 
(D33) 

(D34) 

(D35) 
(D36) 



t+(l/2),i,J+(l/2) 



^i(P)]"7 + (i/2)- (D37) 
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Table 5. Substitutions for Constructing Scalar and Vector van Leer Interpolation 
Functions for Advection in the x 2 Direction from the xi -Advection Definitions 



x\ Advection Function a 


x 2 Advection Function ^ 


Corresponding Substitutions c 


[■*(V)Lj+(l/2) 




[V]*-(l/2)J+(l/2) > M<+(l/2)J-(l/2) 
[H'+(l/2),y+(l/2) > [H'+(l/2),y+(l/2) 

[ 5 l(V)]i-(l/2)J+(l/2) > [ 5 2(V)];+(l/2)J-(l/2) 

[ 5 l(V)]i+(l/2),y+(l/2) > [ 5 2(V / )];+(l/2)J+(l/2) 

[ U l] ( 'J+(l/2) — Ni+(l/2),y 
► My 


L £ ^n (7 uJ i -+(i/2)j+(i/2) 


L^2l02jJ ( - + (i/2)J + (l/2) 


[ a l]ij+(l/2) > [°2]i + (l/2),y 
[ f7 l]/+l,/+(l/2) > [°2]i + (l/2),y+l 

[5l(<Tl)] u+(1/2) ► [52(02)] i+ (i/ 2 ),y 

[^(^Oli+lJ+Cl^) > [^( (7 2)] i+ (l/2)J+l 

[ U l]i,y+(l/2) — » M i+ (l/ 2 ),y 
[ U l]|+l,y+(l/2) > [ U 2]/+(l/2),y+l 
— ► [X 2 ]j 


[</l( <7 l)]i+(l/2)j+(l/2) 




[ ff i]y+(i/2) — > [ a i]u-(i/2) 

[ CT l]/+l,y+(l/2) > [ CT l]/,y+(l/2) 

[ 5 l( CT l)]i,/+(l/2) > [52(<7l)]ij-_(l/2) 

[^l(CFi)]. + 1J+(1/2) ► [5 2 (<Tl)] ii y +( i / 2) 

Mi,y+(l/2) — » M;-(l/2),y 
Ni+l,y+(l/2) > [ U 2] ;+ (i/2),y 

[Xl],- ► [x 2 }j 



a This column contains the van Leer interpolation functions for advection in the x\ direction that have 
been provided in full in the text . 

b This column contains the remaining van Leer interpolation functions not given explicitly in the text. 
These are needed to describe advection in the X2 direction. 

c This column contains the needed substitutions in order to use the x\ interpolation functions 
(eqs. [D9-D1 1] and [D22-D24]) to construct the x 2 functions. 
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Electron-density update: 



At 



^Ji+(l/2)J+(l/2) - ^Ji+(l/2),y+(l/2) [AT 1 ] /+(1/2) 



X (fe] m [^lK/p)]mj + (l/2) - [M-felU^/P)]"?(l/2)) > 

(D38) 

Matter-internal-energy-density update: 

l fi J/+(l/2)J+(l/2) - FJi+(l/2)J+(l/2) [AT!]j +(1/2) 

x (M m [gai] w [^i(^)]"+ij+(i/2) " fe], IM^lS+V)) ' (D39) 

x\ -momentum-density update: 

^1j+(l/2) = N"j+(l/2) " (bli+(l/2) bl] i+ (i/2) [fi (^l/p)]^a/2)J+(l/2) 

" b] i _ (1/2) bl] / _ (1/2 ) ^l(^l/p)]"_ + (T/2)J + (l/2)) , (D40) 

x 2 -momentum-density update: 

[ J2 1"+(V2)J = M^(1/2)J - UyT (W+l [*1 (J2/P)] 

-Mjg3iUM*2/p)]-r)> (° 41 ) 

Radiation-energy-density update: 

At 



L B eJ*+(l/2),i+(l/2),/f(l/2) - L C eJfe+(l/2),i+(l/2),;+(l/2) [ATi] i+(1/2) 



(\g2] i+ i [g3l],-+i [ F l( £ 'e)]i+( a i/2),i+lJ+(l 



/2) 

- [g 2 ] ; foi], [F, (£ e )]^r i/2) ,,. J+(1/2) j • (D42) 



D.3. Update Equations for an x 2 Advection Sweep 

For the x 2 advection sweep we obtain the following set of update equations. 
Density update: 

[p]"-KV2)J+(l/2) = [Pl"+a/2)J+(l/2) - AZ k3l]i+(l/ 2 ) " 

/ b2], +1 [^2(P)]^ /2)J . +1 - [ 8 32]j [H P )Q mj 



[AT 1 ] /+(1/2) [AT 2 ] . +( 



(D43) 



1/2) 



Electron flux: 

Matter-internal-energy-density flux: 

[F 2 ( £ /p)]^ /2)J = 
x\ -momentum-density flux: 

l«+7 



n+7 



'"+(1/2) J 



-i «+7 



^(P)]"+ + (1/2)J 



P / J ;+(l/2) J 



[^('l/p)]U T =[ ( /2(«l)]5 r [^(p)]y r 

^2 -momentum-density flux: 

[f 2 (, 2 /p)]^ /2)j+(1/2) = [M^)t + l m , M m ^(p)]"+ + (i/2)j +( i/2) 

Radiation-energy-density flux: 

-i n+7 



[Wp)C W(1/2)j = (|) 



fc+(l/2),i+(l/2)J 



Electron-density update: 



p+5 
dj+(l/2)J+(l/2) 



l«+7 



^ bl] ; -+(l/ 2 ) ([*l]i+l ~ Wi) 



l ^-+(i/2)J+(i/2) - [AT!]^ [AT 2 ], +(1/2) - 

' (b 2 ] y+1 [F 2 K/p)]^ /2)J+1 - b2],[F 2 (n e /p)]^ /2)J 



Matter-internal-energy update: 



rpin+5 

l fi Ji+(l/2)J+(l/2) 



Az bi],- + (i/2)(^i]i+i-[^i]«) 



jci -momentum-density update: 



M"j~+(l/ 2 ) 



i '+(l/2)J+(l/2) [^l],-+(l/ 2 ) [AT 2 ] i+(1/2) 

x (>32],. +1 [F 2 (£/p)]^ /2) J+1 - fe 2 ]; [F 2 (£/p)]^ /2)! .) 

r nn+y ^ bl]i([^l],-+(l/2)-[^l].--(l/2)) 

[ " lJy+(1/2) [AT!],.[AT 2 ] y+(1/2) 

x ([832] J+l msr/p^l - [g32]j[F 2 ( Sl /p)}^ r ) , 



;t 2 -momentum-density update: 

r in+8 r i't+7 

L J 2Ji+(l/2)J - K- +(1/2)J 



Af[g3l],+(i/2)([*l] m -[*l],-) 

[ATi] i+(1/2) [Ar 2 ] 7 . 

(fe32] y+(1/2 ) [^2(J2/P)]-+J /2)j+(1/2) " [g32]j-_(i /2 ) fc^/p)] -+J /2) 
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Radiation-energy-density update: 

w ,n+8 _ r F in+7 _ * (jf*]m ~ Mi) 

[ fi dfc + (l/2),i+(l/2),;+(l/2) - l fi ej fe+ (l/2), i - + (l/2),;+(l/2) [AT 1 ] /+(1/2) [AT 2 ] ,- +(1/2) 

^ 32 ] 7 - +1 [^(£ £ /p)]Sf 1/2);i . +(1/2)>; . +1 - [g32], [F 2 (£ e /P)C ( 7 1/ 2), ! - + (i/2)j) • (D53) 



D.4. Post-Sweep Update Equations 

After the advected quantities have been updated during an advection sweep, there are several related 
updates that must be performed to render certain related quantities, such as the momentum densities and 
velocities, consistent with one another. The electron fraction Y e should be updated to make it consistent 
with the electron number density n e . The two quantities are related by Y e = n e m,b/p, where is the baryon 
mass. In discretized form, this becomes 

r iff 

\n e 



r yl ffl „ L"gJi+(l/2)J+(l/2) 

^Ji+(l/2)J+(l/2)=^7Tp (D54) 

lPJ/+(l/2)J+(l/2) 



where (0 = a after the first advection sweep and a> = b after the second advection sweep. We must also 
update the velocity field that has changed as a result of momentum being advected. This is accomplished as 
follows 

hh r - 2 ^tj+d/2) 

LPJ/+(l/2)J+(l/2) + LPJi-(l/2)J+(l/2) 

and 

ito Z ^2J»+(1/2)J 

u 2 \i+(i/2) i = 7 ~- (D56) 



[S2]f+(l/2) ([p]/+(l/2)J+(l/2) + [p]i+(l/2)J-(l/2) 

Finally, the temperature and pressure should be related to render them consistent with the updated density, 
electron fraction, and matter internal energy density. The procedure for doing this update is the same as that 
described in appendix I, and we refer the reader there for detail of this procedure. Since the temperature and 
pressure play no role in the advection equations, we perform this update only after the completion of the 
second advection sweep. The equation that must be iteratively solved for the new temperature [T] n+b is: 

E ([ r ]"+(l/2),y+(l/2) . [P]"+a/2)J+(l/2) ' M+*i/2),/f(l/2)) ~ [ £ ]/+a/2)J+(l/2) = °" ( D57) 
Table 6. ADVECTION TIME-SUPERSCRIPT VALUES 

Direction of First Sweep Direction of Second Sweep a j3 y 8 



x 2 



x 2 



a a b 
a b a 
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D.5. The Complete Advective Update 

With the complete prescription for the two advection sweeps now derived and differenced, we are fi- 
nally prepared to assemble the complete sequence of operations of the advection algorithm. This sequence of 
steps is shown in Figure 28 and corresponds to box (a),(b) of Figure 3. The values for the time superscripts 
a, j3, y, and S, which depend on the order of the advection sweeps, are given in Table 6. 



E. Discretization of Viscous Dissipation Equations 

To insure the accurate treatment of shocks, the ZEUS algorithm (Stone & Norman 1992a) includes 
viscous dissipation in the form of a term that mimics the viscous stresses actually present at a real shock 
front. Were such a term omitted from our model, the hydrodynamic equations would cause oscillations to 
develop in the solution. Such oscillations would result from the lack of any means to generate entropy in 
the finite-difference scheme and satisfy the Rankine-Hugoniot jump conditions across the shock. (For more 
detail on the choice and implementation of artificial viscosity schemes, see Bowers & Wilson (1991).) 

In this work we follow Stone & Norman (1992a), but consider only pseudo-tensorial artificial viscosi- 
ties. It has been widely known for sometime (Tscharnuter and Winkler 1979) that tensor artificial viscosities 
offer superior qualities to scalar schemes. In future work, we will consider such prescriptions, but for now 
we restrict ourselves to a scalar approach. Readers interested in tensorial forms for the viscous dissipa- 
tion should consult Stone & Norman (1992a). The ZEUS scheme employs a form of the von Neumann- 
Richtmyer scheme (Noh 1987; Richtmyer & Morton 1967), where the viscous stress in the ith direction is 
given by 

w 2 p(^-\ \£dVi/dxi<0 



Qi=< \dxij (El) 
k if dVi/dxi > 0, 

where w is a characteristic length. With the substitution of l q = w/Axt, where Ac,- is the width of the 
computational zone in the ith direction, we have, upon differencing, 

[0lL- + (l/2)J+(l/2) = l\ [p]" + + (l/2)J + (l/2) {™» " Ny 1 ,0) } 2 (E2) 

[02] /+(1/2)J+(1/2) = l\ [p]£( im j H i/2) {min ([V2]%Ij - [V2\lf ,0) } 2 . (E3) 

By doing this, we have rewritten the characteristic length w in terms of a characteristic zone count. For all 
calculations, we set l q = 2. The quantity l q corresponds to the approximate number of zones over which the 
shock is spread. 



These forms enter into the discretized analogs of equations (39) and (44) exactly as described by equa- 
tions (35)-(37) of Stone & Norman (1992a) and we refer the interested reader there for more detail. We 
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fl(P") [D7] 



Fl(£"/p") [D34] 



£»+»(£») [D39] 



S]»+°(si") [D40] 



F)(s 2 »yp») [D36] 



F^'Vp") [D37] 



£s" + °(£s") [D42] 



F2(p" +n ) [P8] 



F 2 W + <Vp» +n ) D>44] 



£"+*(£"+'■) [D50] 



p»+V + ") [D43] 



F 2 (j 1 "+«/p»+<') [D46] 



Jl"+*(j 1 "+«) [D51] 



E S » + 'W + ") [D53] 



( ENTER ) 




YES """'numbered^ no 



*-| FlW/p") [D33] K 
ntfM«W) [D38] * 



p"+«(P") ff>32] 



F 1 (J 1 "+»/p"+") [D3S] 



J 2 " + °(i2") I 041 ] 



^■W") [D54] 



U^i+n^n+o) [D55] 



U2" +n (52" +0 ) P56] 



;7 e "+*(n c "+») [D49] 



F 2 (£»+a/p/i+o) [D45] 



Y^*Hn^* b ) [D54] 



F2(£ s " +0 /p" + '') [D48] 



F 2 (P") [D8] 



F2(£"/p") [D45] 



£"+»(£») [D50] 



J2 n+<! (J2") P 52 ] 



F 2 (s 1 "/p") [D46] 



F^n/p*) [IMS) 



W ff>53] 



Fi(p"+«) [D7] 



Fl(n e n +<'/p"+n) [D33] 



£»+*(£»+») [D39] 



p»+V + ") [D32] 



F 1 (j 2 " + "/p" + ") [D36] 



S2«+ i (S2"+ a ) [D41] 



F E " +i, (£s"+») [D42] 



F2(n«"/p") [D44] 



nP+W) [D49] 



p»+«(p<>) [D43] 



F 2 (i2"/p") [D47] 



S] n+a( Sl B) [DS1] 



HW) [D54] * 



v l "+"(s l " +a ) [D55] 



u 2 " + "to" +£ [D56] 



M e "+*(n ( />+ ) [D38] 



Fi(£"+<Vp»+<') [D34] 



F 1 (j 1 "+»/p"+«) [D3S] 



Yf+ b {nJ>+ b ) [D54] 



Fi(£ s " +0 /p" +0 ) [D37] 



U2 n+A( J2 n+6) [D56] 



p+6(pn+i£n+4 [D57] 



C ) 



Fig. 28. — A flowchart of the advection sweep algorithm. The blue-shaded sections of the figure show 
the calculations necessary for an advection sweep in the x\ direction. The red-shaded regions show the 
analogous steps for an X2 sweep. The equation numbers corresponding to each step of the algorithm are 
indicated. The initial decision at the top of the algorithm sets up the directional sweeps so that they alternate 
order with timestep. 
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difference equation (39) as 

[ £ ]"+( 1 l/2)J+(l/2) = [ £ ]"+(l/2)J+(l/2) 



A Jr n i ( lvi}i+i,j+(i/2)-l v Aj+(i/2) \ 
-A? |[£2i],- +(1 / 2 )j+(i/2) ^ [xi]j+i _ ^ J 

| [62] l - + (i/ 2 ),j+(i/2) / [^2l,-+(l/2)J+l ~ Ni+(1/2),A 1 E4) 



„ F+i _ r« ^ f [ gl ^+(V2)J+(l/2) ~ [6l]"-(V2) J+(l/2) \ 

u iJ;j+(i/2) - L^iJ^+(i/2) " ftt ftt ~ TTi^+l ri^Tl (b5) 



The components of equation (44) are differenced as 

Ji+(l/2)J+(l/2) ~ [fc?lJi-(l 
M;+(l/2) - [*l]f-(l/2) V [p]"+(l/2)J+(l/2) + [P]h1/2)J+(1/2, 

and 



M+(1/2)J - t U l]"+(l/2)J- 



b] I -+(l/ 2 ) (k] 7 - +( i /2 ) - h] y -_ ( i/ 2 )) 

[62]"+(j/2),j+(l/2) ~ [62]"+(j/2),j-(i/2) \ 
[P]"+(1/2)j+(1/2) + [P]Si/2)j-(1/2) / 



(E6) 



F. Discretization of Gas-Momentum Source Equation 



The solution of equation (43), in substep i, accounts for accelerations due to fluid pressure gradients 
and gravitational forces. Since p remains unchanged during this step of the algorithm, it can emerge from 
the time derivative in equation (43). Thus, we rewrite equation (43) as 



d(v l x l + v 2 x 2 ) 
dt 



= -VP + VO, 

matter P 



(Fl) 



where xi and x 2 are unit vectors in the x\ and x 2 directions. Upon differencing, we obtain expressions for 
updated velocities. In the x\ direction, we have 



[ui\ n+i 



W+(l/2) 



[Ul] 



2At 



n+b 

2At 



'[P 



n+h 



,-in+h 



i+(l/2)J+(l/2) - [ P ]i-(l/2),;+(l/2) 



. [p]?+(l/2) J+(l/2) + [P]"-(1/2)j+(1/2) , 



[ x lJ/+(l/2) ~~ l x Ui-(\/2) 



1/2) V 



At 



[g2 



dx\ 



/2),;+(l/2) 



(l/2)J+(l/2) 



7+1 



(F2) 
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where the last term is the contribution of the time derivative of x 2 . (The time derivative of x\ yields zero 
contribution to this term.) The averaging of the nearest four face-centered values of v 2 in the x\-x?, plane 
gives the value at the need location for calculation of an Vi velocity, viz., at 1,7 + (1/2). 

Similarly, in the x 2 direction, we have 




b] I -+(i/ 2 ) (\x 2 ]j + (y 2 ) ~ [■*2],-_(l/2)) 

(F3) 



p]ti+h Ipin+h 
n i+(l/2),;+(l/2) ~ H,-+(l/2)J-(l/2) 

n+i _i_r n r +1 

(l/2),y+(l/2) + lPJ/+(l/2)J-(l/2) 



Here, we note that for spherical self gravity, the x 2 component of V<I> vanishes. In the case of a general 
gravitational potential, arising from a non-symmetric solution to the Poisson Equation, a non-zero d<t>/dx 2 
would be present. In the x 2 direction, the contribution of the time derivative of xi is absorbed into the 
definition of the velocity. This follows from our definition of u,-, which always has the dimension of length 
over time, even if the dimension of x,-, itself, is one of angle (see §3.2). 

In this algorithm, we assume a 1-D approximation to self-gravity by assuming the mass distribution 
behaves as if it were spherically symmetric. Although not exact, the highly condensed nature of a col- 
lapsed stellar core makes this approximation satisfactory for supernova models. With this assumption of 
a spherically symmetric mass distribution, we obtain a spherically symmetric gravitational potential of the 
form 

,(,,,) __«£(5£>, (F4) 
r 

where M(r,t) is the gravitational mass contained within the radius r at time t. Note that we we have now 
specified spherical-polar coordinates for this calculation. Although, this is the most convenient form for 
problems in spherical geometry, another coordinate system could be used, with the appropriate transforma- 
tion. 

We evaluate M(r,t) as follows. In general, an element of mass is given by 

dM(x,t) = p(x,t)dT 1 dT 2 T 3 , (F5) 

which in spherical geometry is 

dM(r,d,t) =2np(r,d, t)r 2 dr sin Odd. (F6) 

The mass, M(r), is evaluated by performing a radial integral out to the point of interest. When combined 
with an angular integration, we obtain an angularly averaged mass distribution, 

M(r,t) = K f 2np(r,e)r 2 dr&m6dd, (F7) 

02 — 01 JO J 0i 

where the averaging pre-factor in front of the integral allows the computational grid, which runs from 0i to 
02, to be of arbitrary extent — not necessarily reaching either pole. 
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Finally, the differenced version of these expressions can be written. Note, that although the overall 
timestep advancement is from times [t] n to [t] n+l , since a all quantities that go into the gravitational-potential 
calculation are evaluated at time t", the result is also given a time subscript of n: 

[M] ^ = rei. -rei. X X Wlmjnm l^Uw ^W) > ^ 

L Jjmax L Vmin i— /'min 7— imin 

with 

[*]? = -G^. (F9) 

We use / m i n , imax,7min, and y max to indicate lower and upper bounds of jq and X2 dimensions, respectively 
(see Appendix K). Expressions in equations (F8) and (F9) are evaluated once at the beginning of a timestep 
and kept constant until the next timestep. 



G. Discretization of Lagrangean Gas-Energy Equation 

Equation (38) accounts for thermodynamic work done on the fluid internal energy by compression 
and/or expansion. Our discretization is a generalization of scheme of Stone & Norman (1992a), which 
allows an arbitrary equation of state in which the internal energy and pressure are expressed as a function of 
temperature and density. This solution method makes no assumptions about convexity of the EOS. 

We implicitly difference equation (38) as 



E ([ r ]"+( 1 l/2).i+(l/2) ' [p]"+(l/2),y+(l/2) . W"+(l/2)J+(l/2)) ~ [ £ ]"+(l/2)J+(l 



/2) 



+y [ v - y ]"+(i/2)j+(i/2) 

X {p 



jp J+(1/2) , [P] •+( 1 i/2),;+(l/2) > M£(l/2),/+(l/2)) + [ P ]"+(l/2) J+(l/2) } ~ °> ( G1 ) 

where the first and third terms of equation (38) contain the EOS in the form of equations (56) and (57), 
respectively, and where 

f V " V ]"+(l/2 , H/2: TaY~1 ] [ "1 j /+ 1 . /-( l . 2 ! ~ 1X2 j , [X.*, 1 j , .^1 j /. i ■ 2 ) 

P 1 lJi+(l/2) 

4 



JT^- (\g2] i+ i [<?3l] m [Ul]/+i,y+(i/2) ~ [g3l]; N; 

E^^kwi M»<wh - fc^W/»j) • (G2) 



For a realistic equation of state, this equation cannot be solved explicitly, since there is no explicit way of 
expressing matter pressure as a simple function of internal energy. This equation must be solved iteratively 
for the new temperature [T]"^^ j + ^/ 2 y which is unknown and enters parametrically into the new internal 
energy and pressure through the equation of state as described in equations (56) and (57). Since the solution 
of equation (38) is the last step in the algorithm, the new value of the electron fraction and density are 
known prior to the solution of equation (Gl). In practice, we employ a combination of Newton-Raphson 
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and bisection algorithms to accomplish the iterative solution of equation (Gl). The presence of first-order 
phase transitions in complex equations of state, such as those describing matter in core-collapse supernovae 
(Lattimer & Swesty 1991), can create convergence problems for the solution of equation (Gl) Newton- 
Raphson iteration. These same phase transitions also exhibit non-convex behavior, which can also pose 
problems for Riemann solver based methods. When non-convergence occurs during the iterative solution of 
equation (Gl), we make use of the slower converging, but more robust, bisection method for single-variable 
nonlinear equations. 

H. Discretization of Radiation Transport Diffusive/Collision Equation 

H.l. Pair-Coupling of Equations 

In this appendix, we discuss the implicit solution of the sets of diffusion equations for the radiation 
represented by equation (49). These equations are solved as pair-coupled sets described by boxes (c), (e), 
and (g) of Figure 3. By the use of the terminology "pair-coupled," we means that discretized version of 
equation (49) and its antineutrino analog are solved simultaneously for all energy groups for a given type 
of neutrino (electron, muon, or tauon). For example, in the step represented by box (c) of Figure 3, the 
nonlinear equations 

e S e =0 (HI) 

e S e =0 (H2) 

are solved simultaneously over all groups for the electron neutrino energy density e E e and the electron 
antineutrino energy density e E e . There is one such pair of equations for each of the N g energy groups that 
span the radiation spectrum. Note that the diffusion coefficients and the source-terms in equations (HI) and 
(H2) are functions of the neutrino energy and differ between the neutrinos and antineutrinos. The analogous 
pair-coupled sets of equations are solved in the step represented by box (e) of Figure 3 for the muon neutrino- 
antineutrino case and in the step represented by box (g) of Figure 3 for the tauon neutrino-antineutrino case. 

For the sake of compactness of notation in the remainder of this appendix, we use the generic notation 
E e , without the leading superscripts e, jx, or t to refer to the neutrino energy density and E e , without the 
leading superscripts e, /x, or t to refer to the antineutrino density. The finite-difference equations for each 
neutrino species can then be obtained by straightforwardly substituting e E e for E e , e D e for D e , etc. 

H.2. Implicit Finite Differencing of Diffusion Equation 

We discretize equations (HI) and (H2) using the spatial differencing scheme developed by Crank (1980) 
for diffusion equations and using the standard implicit backward-Euler approach for the time evolution 



and 



~ d( e E e 
dt 



' d( e E e 
dt 



. J radiation 



radiation 



-V-(^£ e )-e^rP e :Vv)- 



■ V-r^-e— ( e P e :Vv) 
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scheme. The backward-Euler time discretization scheme has the advantage that method is L-stable for all 
timesteps (Hunsdorfer & Verwer 2003). 



Equation (HI) is implicitly discretized as 

[ £ 'e]^( 1 l/2),i+(l/2)J+(l/2) - [ £ e]^ 



where 



At 

d 

e 



[V • D e V£ e ]£+ (1/ 2) )i+ (i/2) j+d/2) 



<3(P e :Vv)f +1 rc3 



de 



~ [ S e]^+(l/2),i+(l/2)J+)l/2) - ( H3 ) 
J*+(l/2),i+(l/2)J+(l/2) 



[V • AjVE e ]St+(i/2),i+(i/2),y+(i/2) = 
1 



fe] i+ (l/2) [g3l] f +(i/2) [?32]_,-+(i/2) \ [*l]j+(3/2) ~ [ x l 1 ;+( 1/2) 

x ^M,+l [g3l] i+ i [g32];+(i/2) [ £) e(jCl)]t+(i/2),i+lj+(l/2) 

x ^eJ^+(l/2), i -+(3/2)J+(l/2) ~ [ A eJt+(i/2),,-+(l/2),j+(l/2) 
[*l]«+(3/2) ~ [*l]j+(l/2) 

- [g2]j [g3l]j [g32] j+ (i/ 2 ) [ D e(*l)]/t+(l/2),iJ+(l/2) 

x [ £ e]t+(i/2),,-+(i/2)j + (l/2) ~ [ £ e]^+(l/2), i --(l/2).j+(l/2) \ 
[*l]j+(l/2) ~ [ x l]/-(l/2) / 
1 

[*2]/+(3/2) - [ X 2] i '+(l/2) 
/ bl],-+(l/2) [g32jy+i + , 
^ N i+ (l/2) 

J*+(l/2),;+(l/2)J+(3/2) - ^ek+(l/2),i+(l/2)j+(l/2) 



\E, 



[g3l] ; - + (i/ 2 ) [g32]y + , 

b] ; - +(1/2 ) ^WWW 

M+(l/2),/+(l/2)J+(l/2) - [ £ e]S(l/2),i+(l/2),;-(l/2) 



^];+(l/2) " \x2]j. 



(H4) 



(1/2) 



-98- 



, ^(P e :Vv) 

• de 



n+l 



*+(l/2),i+(l/2),/+(l/2) 



L fc J^+(l/2) 

[e]jt+i - M k 



in+1 



X ([ Xe:VV ^+(3/2),i+(l/2)J+(l/2)^r!; (3.2M-+CI 2).J (12. 
- [X e : Vv]^+j 1/2))( . +(1/2)J+(1/2) [£ e ]^(i/2),/ + (l/2)J+(l/2) ) > 



(H5) 



and 



[S. 



in+1 

eJt+(l/2),i+(l/2),y+(l/2) 



L A eJt+(i/2),i+(l/2),;+(l/2) 



1 + - 



17a 



\ 



3 K+(\/2),*+(l/2)J+(l/2) 



V 



[ £ k+(l/2)) 



in+1 



eJt+(l/2),i+(l/2)J+(l/2) L^eJ*+(l/2),i+(l/2),/f(l/2) 



1 + - 



3 [ £ e]fe+(i/ 2 ),i + (l/2),y+(l/2) I [ e ]*+( 



1/2) 



(l £ W(l/2)) 

L [ A£ ]/:+(l/2)[ G ]^+(l/2)i'+(l/2), ( -+(l/2)J+(l/2) ! + " 



[^e]"+(l/2),/+(l/2)J+(l/2) 



([ £ W(i 



/2) 



— C 



1 + 



3 L^eJt+(i/2),i+(l/2)J+(l/2) 



V (J £ W(l/2), 
-1 

L t A£ W(l/2) [ K " S ]i+(l/2)i'+(l/2), 1 -+(l/2)J+(l/2) [ £ 'e]"+(l/2),i+(l/2)J+(l/2) 



iV„-l 



+ c [ £ e]/£(l/2),i+(l/2)j+(l 



/2) 



A/.-1 



L t A£ ]f+(l/2) [ fC " S ]"+(l/2),A:+(l/2). I -+(l/2),j+(l/2) 



f=0 



1 + - 



77a 



V 



(t £ W(l/2)) 



3 [ £ e]"+(l/2), ; -+(l/2)J+(l ( 2i 



(H6) 



In equations (H4), (H5), and (H6) the superscript n + 1 is used to indicate evaluation of the superscripted 
quantity at times n + b, n + d and n+f according to whether the equations for electron, muon, or tauon type 
neutrinos, respectively, are being solved (see Figure 3 for the order of updates). For readability, we have 
omitted displaying all the explicit functional dependences of the microphysics factors (S e , k£, G, K m , K° ut ) 
in the finite-differenced equations. These microphysics factors are also functions of the physical conditions 
of the matter at the spatial point in question (temperature, density, chemical composition, etc.). In addition, 
they are functions of the particle energy or energies in the radiation field. We denote this energy dependence 
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in the discretized quantity through the use of additional subscripts: 

1 n+t 



Mn+t _ 
*+(l/2),i+(l/2),/f(l/2) = 



K a 

[ £ ]«:+(l/2) 



and 



[ ^'1 "+ ( 1/2) ,i+ ( 1 /2) .«+ ( 1/2) ,j+ ( 1 /2) 



i <J J^+(l/2)^+(l/2), ( -+(l/2)J+(l/2) 



i+(l/2)J+(l/2) 
n+t 



C+(l/2)'[ £ ]/t+(l/2)) 
G ([ £ ]c+(l/2) ) t £ ]/t+(l/2) 



«+(l/2)J+(l/2) 
,-+(l/2)J+(l/2) ' 



Equation (H3) is fully implicit in terms of the unknowns [E t 



m+l 

J*+(l/2),,+(l/2)J+(l/2) 



(H7) 
(H8) 

(H9) 
and 



[^e]fc+(i/2) (+(i/2) y+(i/2) ^ Ut some °f tne coefficients arising from the microphysics are evaluated at the 
beginning of the substep. Following the practice of Turner & Stone (2001), we evaluate the diffusion co- 
efficient and Eddington tensor based on information at the beginning of the substep. The reason for this is 
that the absolute value function that appears in the evaluation of the Knudsen number would introduce a 
discontinuity in the Jacobian of the nonlinear system if evaluated implicitly. Such discontinuities can often 
produce failures in Newton-type iterative procedures for nonlinear systems. By evaluating the diffusion co- 
efficient at the beginning of the substep, we avoid this difficulty altogether. This approximation is accurate 
as long as the change in the radiation energy density is small over a timestep. For most astrophysical heat- 
ing/cooling problems, this is indeed the case. In equation (H3), all quantities that are temperature-, density-, 
or 7 e -dependent are also evaluated at the beginning of the substep. Doing so avoids the difficult problem 
of simultaneously solving the radiative diffusion equations while also solving the Lagrangean portion of 
the gas-energy equation. All quantities, such as the opacities, emissivities, etc., are evaluated by using the 
values of the temperature, density, and electron fraction that are known at the beginning of the substep. If 
the exchange of energy, or lepton number, that takes place during a timestep is small, this approximation is 
accurate. The validity of these approximations are tested in the verification problems that we consider in 
§4- 

The antineutrino counterpart of equation (H3) is given by 



^eE(i/2),i+(l/2)J+(l/2) - [£e]*+(l/2),,+(l/2)J+(l/2) 



At 



[V • fieV£ e E( 1 /2) )i+ (i/2),y+(i/2) 



d (P e : Vv) 
de 



n+\ 



[^ e ]/t+(l/2),i 



*+(l/2),i+(l/2),/+(l/2) 



*+(l/2),i+(l/2)J+)l/2) 



= 



(H10) 
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[v • o e V£ e ]^ (1/2))i . +(1/2) = 



i i 



x 



[S2]/+(l/2) [^3l]/+(l/ 2 ) [§32] 7 - + (i/ 2 ) \ [ x l] i+ (3/2) ~ i x l]i+(l/2) 

[g2] i+ i [g3l] i+ i k32]y+(i/ 2 ) [AK*l)]*+(l/2),mj+(l/2) 

x [^e]t+(i/ 2 ),,-+(3/2) J+(l/2) ~ [^e]^+(l/2), i -+(l/2)J+(l/2) 
[*l]i+(3/2) - [*l]i+(l/2) 

- [g2] ; [g3l]; [g32] j+{l/2 ) [De{xi)} n k + {l/2) jj +{l/2 ) X 
^eJfc+(l/2),i+(l/2)J+(l/2) - L£ek+(l/2),i-(l/2)J+(l/2) \ 



[ X l]i+(l/2) - [*l]f- 



(1/2) 



X 



[ X 2] 7+ (3/2) - [ X 2] i '+(l/2) 
L?3l]; + (i/ 2 ) [gn] j+l n+t 

M i+(1/2) PtW W)WWi 

x [^e]t+(i/ 2 ),,-+(l/2) J+(3/2) ~ [^e]^+(l/2), i -+(l/2).j+(l/2) 
[ x 2] i -+( 3 /2) - [ x 2]y+(i/ 2 ) 

[g3l],-+(l/2)[g32]y + , 

x [^e]fc+(l/2),i+(l/2) J+(l/2) ~ [^e]^+(l/2), i -+(l/2).j-(l/2) 
[ x 2] i -+(i/ 2 ) - [ x 2] 7 -_(i/ 2 ) 



5(P e :Vv) 

e 



77+1 



[ £ W( 



1/2) 



[gl _ [gl 

*+(l/2),i+(l/2),/f(l/2) L J * +1 L Sk 



ds 

[X e : Vv]^J 3/2) ; . +(1/2)j . +(1/2) [£e]^+( 1 3/2),i+(l/2)J+(l/2) 



[X e : Vv]^J 1/2))( . +(1/2)j . +(1/2) [£e]^(i/2),i+(l/2)J+(l/2) 
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and 



[S. 



n+t 

£Jfc+(l/2),/+(l/2)J+(l/2) 



1 + 



3 [^e]K(i/2),i+(l/2),j+(l/2) 



S 1"+! 
fc+( 

+ c [^]*+(l/2),i+(l/2),j+(l/2) [£eJ£+(l/2),i+(l/2)j+(l/2) 



e U+(l/2),i+(l/2).i+(l/2) 
"P. 



\ 



([ £ W(l/2)) 



in+1 



~3 ^eJfe+(l/2),i+(l/2),y+(l/2) 



([ e ]*+(l/2)) 



[ £ W(l/2) 



X L I A£ W(l/2) [ G ]"^(l/2)^+(l/2), ( -+(l/2),;+(l/2) 
(=0 



1 + - 



T7a 



3 K+(l/2),i+(l/2)J+(l/2) 



([ £ W(l/2)) 



-C 1 + 



([ £ W(l/2)) 



[^e]jfc+(l/2), 



i+(l/2)J+(l/2) 



x £ [Ae 



(=0 



+ c [^e]^(l/2),i+(l/2),y+(l/2) 



P. 
/ 



p ,n+l 



^+(l/2)[ ,cS ]i+(l/2)i+(l/2), ( -+(l/2)J+(l/2)^ e j£ + (l/2),i+(l/2)J+(l/2) 



X L t A£ W(l/2) [ S '1"+(l/2)^+(l/2),;+(l/2)J+(l/2) 
f=0 



1 + - 



([ £ W(l/2)) 



3 L^eJ/+(i/2),i+(l/2) J+(l/2) 



(H13) 



The finite-differenced antineutrino equations (H10)-(H13) are almost analogous to the finite-differenced 
neutrino equations (H3)-(H6). The one exception to this one-to-one correspondence is in the pair-production 



terms where the same discretized function [G] 



n+t 



<+(l/2),t+(l/2),i+(l/2) ) j+(l/2) 



appears in both equation (H6) and 



equation (H13). In the latter, the summation is over the first index of the function, while in the former, the 
summation is over the second index. Because both E e and E e appear as factors in this set of terms, the 
process of pair-production is responsible for coupling the neutrino diffusion equations to the antineutrino 
equations. This requires that the two sets of equations be solved simultaneously for each of the three neutrino 
flavors: electron, muon, and tauon. We discuss the method of solution of these implicitly differenced 
equations in §H.6. 



H.3. Implementation of Boundary Conditions 



Boundary conditions are specified for equations (HI) and (H2) by altering the differencing of the V • F e 
and P e : Vv terms. In this subsection, we only present boundary conditions corresponding to equation (HI), 
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with the understanding that similar boundary conditions for equation (H2) can be trivially obtained by 
substitution of E e for E e etc. 

We consider three types of boundary conditions in the x\ direction: Dirichlet boundary conditions (E e 
is specified), zero-flux boundary conditions, and free-streaming boundary conditions. The first of these 
conditions is handled by equation (H4) with no modifications necessary. Zero-flux boundary conditions are 
usually applied at the inner or left-edge of the grid. To implement zero-flux boundary conditions on a given 
edge, the corresponding DWE terms in equation (H4) for zones at that edge are set to zero. In the case 
of free-streaming boundary conditions at the outer (or right) edge of the grid, the corresponding DWE is 
replaced by cE e in an upwinded sense. Equation (H4) then replaced in the outermost zone by 

[V-D e V£ e ]^ (1/2))i . +(1/2)j . +(1/2) = 



[S2]/+(l/2) L?3l],- + (i/2) [g32]j+(i/ 2 ) [ [*l]/+(3/2) ~ [ x l],-+(l/2) 

X {[gl\i+\ L?3l] m k32]y + (i /2 )C[£ e ]^ ( 1 1/2) ,,- +(1/2)J+(1/2) 

- \g2\i [gZ\\i [g32]y+(i/ 2 ) [ D e (*1 )]fc+(l/2),iJ+(l/2) 

x [ £ e]t+(i/ 2 ),,- + (i/2)j+(i/2) ~ [ E eTkt(i/2)J-(l/2),j+(i/2) \ 
\xi\i+{l/2) ~ [ x l]i-(l/2) / 
1 

[•*2]j + (3/2) - i x 2]j+(l/2) 

( tg3l],- + (i/ 2 ) [g32]y+l , , , ]n +t 
x ^ b];+(i/2) [Ds(x 2 )] k+{1/2)M1/2)J+1 

x [ A eJt+(l/2),i+(l/2)J+(3/2) ~ ^ek+(l/2),,-+(l/2)j+(l/2) 
[ x 2] i -+( 3 /2) - [ x 2]y+(i/ 2 ) 

[g3l] i+{1/2) [g32]j +t 

[ L, e{X2)\ k+{l/2)Ml/2 )j 



[g2\i+(l/2) 

PeJ*+(l/2),,-+(l/2)J+(l/2) - [ fi eJfc+(l/2),i+(l/2),;-(l/2) 



[ X 2] i '+(l/2) - 



(H14) 



In this equation, we have made the assumption that radiation is flowing off the grid at the outer edge in order 
to determine the upwind direction. This boundary condition must be used with care and only in cases where 
one is sure of the accuracy of this assumption. 

In the x 2 direction, we consider three types of boundary conditions: Dirichlet, zero-flux, and periodic. 
The first of these is again trivial and requires no modification of equation (H4). The second, zero-flux 
conditions, is treated as we previously described, by setting the appropriate DVE terms to zero in equation 
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(H4) for the zones at the corresponding edge of the mesh. The third case, periodic boundary conditions, is 
treated by using the ability of MPI to establish periodic process topologies. Using MPI send and receive 
subroutines to exchange ghost zones, together with a process topology that is periodic in the x 2 direction, 
permits equation (H4) to be used to apply periodic boundary conditions without modification. 

Spectral boundary conditions for the P e : Vv term in the e direction are handled similarly to the pre- 
viously described zero-flux spatial boundary conditions. The appropriate terms of P e : Vv in equation (H5) 
are set to zero at the lower and upper edges of the energy grid. This ensures that no energy flows out of the 
spectrum at either end. 



H.4. Evaluation of Diffusion Coefficient and Eddington-Tensor Terms 



In the verification tests described in this paper, we rely on the Levermore & Pomraning (1981) pre- 
scription for flux-limiting, which relates the diffusion coefficient and the Eddington tensor to the Knudsen 
number for the radiation flow. However, the modification of our scheme to utilize other flux-limiters is 
trivial so long as the diffusion coefficient and the Eddington tensor can be cast into a form that relies on the 
Knudsen number. 

In what follows, we will generically describe the evaluation of Knudsen numbers, diffusion coeffi- 
cients, and Eddington tensors for neutrinos. The analogous equations for antineutrinos can be obtained by 
straightforward substitution of E e for E e , k a for k°, etc., and is not presented here for reasons of brevity. 

In our calculation of the diffusion coefficients and Eddington tensor, we emulate Turner & Stone (2001) 
and make use of separate Knudsen numbers R e l and R e2 to describe the radiation flow in each of the two 
orthogonal coordinate directions. Since Knudsen numbers are based on components of the gradient of 
radiation energy density, which is defined at cell-centers, the gradients, and hence Knudsen numbers, are 
naturally evaluated at cell faces — precisely where they are needed. First, for /? £l ,we have 



[R 



el. 



*+(i/2),y+(i/2) 



|V£ e -xi| 

KjE e 



n+t 



*+(l/2),i,/f(l/2) 

expanding the right-hand-side we can write R £ l as 



1 



dE F 



n+t 



*+(i/2),y+(i/2) 



(H15) 



m+t 



eU*+(l/2),y+(l/2) 



c e Jt+(l/2),i+(l/2)J+(l/2) 



+ [*e]t+(l/2),i-(l/2),/f(l/2). 



\rp ]tl+b , r-rp m+b I 

. ^£J*+(l/2),,+(l/2)J+(l/2) + l fi eJfc+(l/2),Hl/2)J+(l/2) / 
^eJ*+(l/2)+ ; -+(l/2)J+(l/2) - L A eJt+(i/2),i-(l/2)j+(l/2) 



[ x i];+(i/2) ~~ t x i]i-(i/2) 



(H16) 



In equation (H16), we have arithmetically averaged the opacities and radiation-energy densities in adjacent 
cells to obtain values at the zone interface in the x\ direction. The most recent values of the radiation-energy 
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densities, obtained at the end of the advection substep, are used in equation (H16). The values for the 
opacities are obtained at the beginning of the substep at point n + t, where t takes on the values b, d, or / for 
electron, muon, and tauon type neutrinos respectively. In the X2 direction, an analogous expression can be 
written 



[*e2]*+(l/2),i+(l/2)j 



[J 1 1 fl~T~t r r j'~\ H~\~t 

K e W(l/2),i+(l/2)J+(l/2) + l K e Jt+(l/2),i+(l/2),y-(l/2) 



n+b , \rp yi+b 

eJt+(l/2),i+(l/2),j+(l/2) + l fi eJfe + (l/2),i+(l/2),;-(l/2) , 



n+b 



l/2),/+(l/2)J+(l/2) 



i E e]kt(i/2)J+(l/2),j-(i/2) 



[82 



*+(l/2) 



y+(i/2) 



(1/2) 



(H17) 



To obtain the diffusion coefficients in each direction, D e (x\) and D e {x2), we note that D e , in the 
Levermore-Pomraning formalism, is given by equation (11) as 



D E = 



cX e (R e ) 



KL 



(HI 8) 



where A e is as given in equation (12) and fcj is as given in equation (10). Many other flux-limiting pre- 
scriptions can be cast into this form and could easily be accommodated in this algorithm by replacing 
equation (12) with some other function (see Janka (1991) for a description of how numerous flux-limiters 
fit this form). Thus, we have: 



t+(l/2),i,7+(l/2) 



cX(R e 1 +t ) 



n+t 

J*+(i/2),y+(i/2) 

n+t 



l/2),y+(l/2); 



[*«]*+( l/2),i+(l/2),;+(l/2) + [ K l\k+\ l/2),i-(l/2),;+(l/2) 



(H19) 
(H20) 



and 



[De(x2)r' 



t+(l/2),i+(l/2),; 



£2J 



n+r 

*+(l/2),i+(l/2),; 



2cA([7? e2 ]^ + ( 1/ 2) ;i - + ( 1 /2)j) 



[ ,C e 7 ']i+(l/2),i+(l/2).i+(l/2) + [ K e]fc+(l/2),i+(l/2),;-(l/2) 



(H21) 
(H22) 



The evaluation of the Eddington tensor components in terms of the Knudsen number is similarly 
straightforward. Equations (16) and (17) describe the Levermore-Pomraning prescription for the Eddintgton 
tensor. The algorithm we present is trivially modifiable to handle other closures replacing these equations 
with other alternative functions describing the tensor components in terms of the Knudsen number. 
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We now focus on the components of the Eddington tensor needed to evaluate the tensorial double- 
contraction terms of equation (H5). This double-contraction can be expressed as 



X e :Vv = {XJ n {Vv} n + {XJ 12 ({Vv} 12 + {Vv} 21 ) 
+ {X e } 22 {Vv} 22 + {X e } 33 {Vv} 33 . 



(H23) 



We have dropped terms that are identically zero in a two-dimensional formulation, namely, those propor- 
tional to {X e } 23 {Vv} 23 and {X e } 32 {Vv} 32 . 

The components of Vv are evaluated using the expressions derived in Appendix H.5. The required 
components of X e are evaluated as follows. Using equation (15), we have 

n+t 



[{X e }n 
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,(H24) 
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,(H25) 



(H26) 



(H27) 



J t+(l/2),i+(l/2) J+(l/2) 

The scalar Eddington factor # e is evaluated using equation (16), which requires Knudsen numbers, which 
we have evaluated at cell faces. To perform the evaluation of % e , which is required at cell centers, we average 
over the surrounding (directionally decomposed) face-centered values, viz. , 

K^(l/2),i+(l/2),y+(l/2) = 7 



[Xe (1 +^e(^el) 2 )] (t +( 1 / 2);i - + ij + ( 1/2) 
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+ [A e (l+A e (/? e2 ) 2 )]^ (1 



fe+(l/2),i+(l/2) J J ' 



(H28) 



where we use equation (12) in our evaluation of the A e 's. 



The factors in equations (H24-H27) containing VE E and dE e /dxj arise from the definition of the radi- 
ation flux direction n. They are evaluated at cell centers by centered differences: 
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(H29) 
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and 
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(H30) 



(H31) 



With the preceding derivations, we now have complete information regarding the differencing of the 
velocity-dependent terms. 



H.5. Velocity Gradients 



The evaluation of equation (H23), described in the previous subsection, requires the gradient of the 
velocity field, Vv. Since Vv is a second-rank tensor, we can use equation (126) in Appendix A of Stone & 
Norman (1992a) to help with the evaluation. Upon inspecting that equation, it is easy to extract the set of 
non-zero elements of Vv and obtain 



{Vv} 



33 



gl dx 2 g 2 dx\ ' 

Ul dg31 U 2 dg 32 



531 dxi g 2 g32 dx 2 



and 



{Vv} 21 



1 dV\ V 2 dg 2 
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When differenced, these equations yield the following expressions for the three diagonal elements, 

N;+lJ+(l/2) - [Wl],-J+(l/2) 
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(H32) 
(H33) 
(H34) 
(H35) 

(H36) 
(H37) 
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and 
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(H39) 



2 [g2],+ (l/2) [g32]y + (i/2) 

The evaluation of equation (H23) requires that components of Vv be evaluated at cell centers. Although this 
is a natural location for evaluation of the diagonal components of Vv, this evaluation location does not fall 
naturally for {Vv}i2 and {Vv} 2 i- Because of cross derivatives, these off-diagonal terms are more naturally 
evaluated at the vertices of the integer mesh. To obtain the values at the needed locations, we evaluate these 
expressions using a four- way arithmetic average of the surrounding vertex values, i.e., 

[{Vv} 12 ] i . +(1/2)J+(1/2) = \ ([{Vv} 12 ]„ . + [{Vv} 12 ], . +1 + [{Vv} 12 ] . +1J + [{Vv} 12 ] mj+1 ) , (H40) 
which, upon use of equation (H35), gives 
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Following the same procedure for {Vv} 21 , and using equation (H36), we arrive at 
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H.6. Iterative Solution of the Implicit Finite-Difference Equations 

The implicit finite-differencing scheme for the numerical solutions of equations (HI) and (H2) is fully 
described in the previous three subsections of this appendix. This finite-differencing results in a coupled 
set of nonlinear algebraic equations that must be solved for the complete set of radiation energy densities 
for each flavor of radiation. In this section, we briefly describe the iterative method employed to find the 
solution of this set of nonlinear equations. 

The complete set of implicitly discretized equations that must be solved for each radiation flavor is 
given by equations (H3) and (H10). The unknowns in this set of equations are [£ e ]^,j , 2 > i+ n/ 2 ) y+(i/2) an< ^ 

[■^£]jt+(l/2) ('+(1/2) y+(l/2) ^ 1) • • • >^Ys>> i 'min) • • • j'max> j Jminj • • • jj'max- W(3 USe i mm , 2 ma x > Jmin > <ind 

Jmax to indicate lower and upper bounds of x\ and x 2 dimensions, respectively (see Appendix K). 

Nonlinearities in this set of equations arise in two forms. First, the pair-production terms in equations 
(H6) and (HI 3) give rise to a bilinear coupling between neutrinos and antineutrinos if the pair-production 
kernel G is non-zero. This is what necessitates the simultaneous solution of equations (H3) and (H10). The 
second set of nonlinearities are quadratic in form and arise from the Pauli-blocking factors that appear in 
the scattering integrals in equations (H6) and (H13). In the case where a = the nonlinearities disappear 
and the set of equations decouples into two separate sets of linear equations: one set describing the energy 
densities of the radiation particles and one set describing the energy densities of the radiation antiparticles. 

In the nonlinear case where a and G are non-zero, this set of equations could be solved by various 
means. Because of the sparseness of the nonlinear system, we have used a Newton-BiCGSTAB variant 
of Newton-Krylov iteration (Kelley 1995). This iterative algorithm consists of an outer Newton iteration 
(Kelley 2003) loop surrounding a Krylov subspace iteration (Saad 2003) to solve the linearized Jacobian 
system. The overall flow of the algorithm is depicted in Figure 29. The BiCGSTAB iteration (Kelley 
1995; Saad 2003) is implemented exactly as described in Barrett, et al. (1994). The iterative convergence 
of this method is hastened, or enabled, through the use of preconditioners (Chen 2005). The particular 
preconditioning strategy that we employ is a 2-D extension of the sparse approximate inverse approach used 
for 1-D multigroup flux-limited diffusion problems (Swesty, Smolarski, & Saylor 2004). 

This iterative Newton-BiCGSTAB approach approach has several distinct advantages. First, the matrix 
corresponding to the Jacobian need not be computed and stored; instead, the Jacobian can be applied in 
operator form. Second, the knowledge of effective preconditioners from the linear flux-limited diffusion 
equation can be exploited to accelerate iterative convergence. Third, the method, including preconditioning, 
is easily parallelizable. 

One minor deviation from standard Newton-Krylov iteration that we employ is the enforcement of 
positivity during the nonlinear iteration step. In Figure 29, we indicate this by the red box prior to the 
solution update. During this step, the Newton step for a particular variable is limited if it would result in a 
negative value for the energy density. 
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Fig. 29. — A flowchart of the Newton-Krylov iteration algorithm used to solve the implicitly discretized 
nonlinear radiation diffusion equations. In this flow chart, for simplicity, we have omitted error conditions, 
such as non-convergence of either the Newton or the Krylov loops. Such conditions are always considered 
as fatal, and the calculation terminates. 
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The stopping criterion for the Newton iteration is 

max | 1 1 | < £ Newt , (H43) 

where the solution tolerance, £Newt, is chosen to be some value between 10~ 6 and 10~ 8 . In equation (H43), 
is the £th estimate of the unknown vector is given by 

X T = (e»+\. . . ,E»+\E? +1 ,. . . \) . (H44) 

The stopping criterion for the Krylov iteration (BiCGSTAB) is based on the norm of the linear system 
residual r, and is given by 

M < £Krylov|b|, ( H45 ) 

where b is the right-hand-side of the linearized system of equations on each Newton iteration. We typically 
choose £Kryiov = 10~ 2 £Newt, based on extensive numerical trials, to yield an inexact-Newton method. 



I. Discretization of Energy and Lepton Exchange Equations 

In the radiation-transport equations in Appendix H, collisional processes that change lepton number 
and/or energy of the components of the radiation field are included in the terms S e and S e in equations (HI) 
and (H2) respectively. Since these terms depend on the radiation energy densities E e and E e , once equa- 
tions (HI) and (H2) have been solved, the amount of energy and lepton (if any) exchange is thus determined. 
Once the values of S e and S e are known, equations (35) and (40) can then be solved to accomplish steps d, 
f, and h of Figure 3. 

The right-hand side of equation (35) gives the total amount of lepton number exchange with a given 
flavor of radiation in terms of an integral over the spectrum: 

N = -l f — ^ — -Us. (II) 

The leading superscript e on quantities in equation (II) is used to indicate radiation of the electron neutrino 
flavor (muon and tauon type neutrinos do not contribute to lepton number exchange since they are only 
produced in particle-antiparticle pairs). The first term in the integrand of equation (II) gives the lepton 
number exchange between matter and radiation particles while the second term gives the lepton number 
exchange with antiparticles (if they are present). The minus sign in the second term reflects the negative 
quantum number assigned to the antiparticles. 

In our discretization of equation (35), we replace the integral with a midpoint-rule summation and the 
time derivative with a forward difference to get 
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which is the discretization of equation (35). Once the new value of the electron number density is known, 
the new electron fraction can be computed trivially as 
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+(l/2),;+(l/2) 
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The solution of the energy exchange equation (eq. [40]) is similar to the solution of equation (35), but 
involves several additional steps. Although lepton number exchange only takes place for neutrinos of the 
electron flavor, the exchange of energy takes place for all flavors. The right-hand-side of equation (40) is 
given by 



S = 



(14) 



where the leading superscript I is used to indicate the flavor of the radiation-species (in the case of neutrinos 
I equals e, jj, or t). Equation (40) thus becomes 

~dE~ 

J collision 



dt 



(15) 



The concept of operator splitting can be applied to equation (15), so we can solve separate equations of the 
form 

' d 4)j =- f(%+ (: S e )de (16) 

at / U collision-** J V 7 

for each radiation species. These operator split sub-equations are each solved immediately after the solution 
of the radiation diffusion equation for the corresponding radiation flavor. In the case of electron type neu- 
trinos in substep d (box d of Figure 3), equation (16) is solved after equation (12) is solved. Equation (16) is 
solved for muon neutrinos in substep / and for tauon neutrinos in substep h. 

The discretization of equation (16) proceeds almost identically to the discretization of equation (II). In 
substep d the discretized equation solved to account for energy exchange with electron neutrinos is 
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In substep / the discretized equation solved to account for energy exchange with muon neutrinos is 
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In substep h the discretized equation solved to account for energy exchange with tauon neutrinos is 
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In each of substeps d, f, and h, it is necessary, after the new value of E is obtained, to use the equation of 
state to obtain new values of the temperature and the pressure corresponding to the updated internal energy 
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density and, in the case of substep d, the new updated electron fraction. This is performed by inverting the 
EOS to find the temperature corresponding to the new energy density via a combination of Newton-Raphson 
and bisection methods for nonlinear equations. 

In substep d the new temperature is found by iteratively solving 

E ([ r ]"+(1/2)J+(l/2) . iP]'i+(l/2)J+(l/2) > [ Ye ]'i+(l/2),j+(l/2)) - [ £ ]"+(1/2)J+(l/2) = ( I10 ) 

for the new temperature, based on the updated values of E and Y e . In the case of substep / the equation 
solved for the new temperature is 

E ([ r ]"+(l/2),;+(l/2) ' [P]"+( 1 l/2)J+(l/2) ' [ 7e ]"+(l/2),;+(l/2)) - [ £ ]"+(l/2)J+(l/2) = °' t 111 ) 

while in substep h the equation solved is 

E ([ r ]"+(l/2)j+(l/2) . [p]"+(l/2),;+(l/2) . [ Ye ]'i+(l/2),j+(l/2)) ~ [ E ]"Z(l/2),j+(l/2) = °- ( I12 ) 

In each case, once the new temperature is known, the pressure is determined by the equation of state. The 
computational cost of the equation of state computations needed to solve iteratively equations (I10)-(I12) is 
negligible compared to the computational effort involved in solving the nonlinear systems corresponding to 
the radiation diffusion equations. 



J. Discretization of the Gas-Momentum Radiation Equation 

In this appendix we describe the solution of equation (45). Following the update to gas momentum 
from the matter contribution, as described in Appendix F, we next update the gas momentum to account 
for radiation-matter interaction. This update corresponds to the solution of the [<3(pv)/d?] radiation term in 
equation (41), 

= -V-P rad . (Jl) 

radiation 

For radiation transport algorithms that have a non-zero radiation-matter momentum exchange, this step of 
the algorithm can also be used to solve for the [d{p\)/dt] conhion term in equation (41). In such a case, the 
right-hand side of equation (Jl) will have an additional term P. As previously stated in §2, we assume P = 0. 
(However, see §2.4 for a discussion of radiation-transport algorithms and non-zero value of P.) 

To evaluate V- P ra d, we use Stone & Norman (1992a), equations (130)-(132). Note that several terms 
of these equations in Stone & Norman (1992a) have sign errors. Specifically, the dh\/dx\ term in their 
equation (130) and the dfi2/dx2 term in their equation (131) both have the incorrect signs. Fortunately, 
these errors are irrelevant since the choice of the hi and hi for the coordinate systems under consideration 
makes these terms vanish. In writing the following expressions, we have omitted terms in the divergence that 
are always zero in our implementation. These include all terms in d / dxj, , terms proportional to zero-valued 
tensor elements P13, P23, terms proportional to derivatives of the constant function h\, and terms containing 



d(v 1 x l + v 2 x 2 ) 
dt 
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(J2) 



dg^\/dx2 and dgsijdxx, which are zero by definition. The remaining non-zero terms give 

{V-P rad } (1) = { ^— (g2g3l832 {Prad} n ) + 4-(g3lg32{P r ad} l2 ) 

y ' g2g3\g32 L OX\ dX2 

{Prad} 22 d §2 _ {P r ad} 3 3 d (g3lg32) 
g2 dx\ g3lg32 dx\ 

{v-p rad } (2) = ^{^(^ 32 {P ra 4 21 ) + ^(^f {P™4 22 

{Prod} 33 d (g3lg32) | ({ ? rad] n + { ) dg 2 ^ 

g2g3lg32 dx 2 g2 dxi 

{V-P rad } (3) =0. (J4) 

The quantity P rad represents the radiation pressure tensor containing contributions from the complete 
spectrum of radiation. Hence, 

Prad = L^P e (+*P £ )rfe, (J5) 

where the sum over / represents contributions from all species of radiation. The added term in parentheses 
indicates the possible presence of a distinct antiparticle that is also evolved. Remembering the definition of 
P ra d (see eq. [14]), this expands in terms of the Eddington tensor to give 

P -d = L_// X A (+%Ee) de, (J6) 

where, once again, the possible presence of an antiparticle is indicated in parentheses. 

Expressing equation (Jl) in differenced form (remembering that the density p is constant over this 
step in the operator splitting and emerges from the time derivative), we have the following, where here, the 
possible antiparticle has been suppressed to make the equations readable, and the radiation-energy spectra 
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have been discretized using index k. This update corresponds to substep j (see §3.4): 

2At 

■ - ,,. " ■ ! 

J ;j+(i/2) 



7+ (1/2) i^n+l 



III 



[P]"+(l/2)J+(l/2) + [Pli-(l/2)j+(l/2) 
1 



* Vtelibl]; ([jCl]/+(l/2) - [^l]/-(l/2)) 
x{H ; - + ( 1/2 ) bl] i+( i/ 2 ) f p{Xe} n 

'{X e } n 



t+(l/2),i+(l/2)J+(l/2) 



[ £ eW(l/2),;+(l/2),y+(l/2) 



- M,--(l/2) [^3l]f_(l/ 2 ) 
1 



*+(l/2),i-(l/2),/f(l/2) 



[ £ e]^+(l/2),,--(l/2)J+(l/2) 
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{X e } 12 
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+ 
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£ { X ell2 



*+(l/2),i+(l/2),y+(3/2) 



*+(l/2),i-(l/2),/f(3/2) 
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[ £ eW(l/2),i+(l/2)J+(l/2) 



[ £ e]jt+(l/2),i- 



(l/2)J+(l/2) 



b 2 ] y (f{x e } 12 

{ X e}l2 



+ 



+ 



+ 



! {x e } 12 



t+(l/2),i+(l/2)J+(l/2) 
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+ 



2 L?3l]i L Ji 



+ 



{X e } 22 

1 {Xe} 33 

{X e } 33 



*+(l/2),i-(l/2),y+(l/2) 



[ £ 4+(l/2),i-(V2)J+(l/2) 



fe+(l/2),i+(l/2)J+(l/2) 



t+(l/2),i-(l/2),y+(l/2) 
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+(l/2)j+(l/2) 



[ £ e]^+(l/2),--(l/2)J+(l/2) f ( J7 ) 
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2M 



[ U 2],- + (i/ 2 )j n+l , f , 
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In this two-dimensional formulation, U3 = at all times. The summations run over energy index k and 
species index i. These expressions are evaluated over the entire spatial domain, i.e., for each index i and j. 
To include the effects of antiparticles in equations (J7) and (J8), the following substitution should be made 
in equations (J7) and (J8) for each instance of the expression to the left of the arrow: 



{X e }, 



fe+(l/2),i+(l/2)J+(l/2) 

'{X £ }„ 



? 4+(l/2),i+(l/2)J+(l/2) 



+ 



1 y h 



*+(l/2),i+(l/2)J+(l/2) 



t+(l/2),i+(l/2)J+(l/2) 



-eJt+(i/2),i+(l/2)J+(l/2) 



s]*+(l/2), 



i+(l/2)J+(l/2) • 



(J9) 



K. Implementation of Explicit Boundary Conditions 

In this appendix we discuss the formulation and implementation of boundary conditions for the equa- 
tions that are solved explicitly. Boundary conditions for the implicitly solved radiation-diffusion equations 
are discussed in Appendix H. 

To permit solution of a wide range of problems, we have implemented multiple options for boundary 
conditions into V2D. In this appendix, we describe four such options: (i) flat, (ii) periodic, (iii) reflecting, 
and (iv) flux conserving. V2D makes no requirement that the same option be chosen in each direction or at 
each of the two boundaries in a given direction — in many problems, it is common for them to be different. 

As shown in Figures 30 and 31, V2D is structured to use boundaries that are centered on the integer 
mesh, i.e., boundary interfaces lie on cell faces, not cell centers. In this appendix, for each option listed 
above, we present the necessary boundary-value assignments for both scalar and vector quantities. Since our 
differencing scheme uses values no farther from the point of interest than that of next-to-nearest neighbor, 
we need carry no more than two boundary cells. 

In what follows, we describe the application of boundary conditions only in the x\ direction. Appli- 
cation of the same sets of conditions in the X2 direction is straightforward, as long as one remembers that 
components of vectors normal and tangential to the boundary interface are swapped. 

1. Flat Boundary Conditions. To implement flat boundary conditions in the x\ direction, we apply the 
following procedure: First, for any scalar quantity y, at the left-hand boundary, we set the values at the two 
left-hand boundary zones to the value at the first physical zone: 

[V] !min -(3/2)J+(l/2) = [H mm -(l/2)J+(l/2) = M !mi „ + (l/2)../+(l/2) ■ ( K1 ) 

This applies to all values of the %2 index j, which is also the case for all such equations in this appendix. 
Similarly, to implement flat conditions at the right-hand boundary, we set the scalar values of the two right- 
hand boundary zones to have the value of the last physical zone: 

M imax +(5/2)J+(l/2) = Mw+(3/2)J+(l/2) = [V] imax +(l/2)J+(l/2) ■ ( K2 ) 
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•/min + ^ 
Jnuxi + ^ 



'min ~ ^ 'min ~ ^ 'min 'min + ^ 'min + ^ 



Fig. 30. — The bottom left-hand comer of the staggered mesh. The shaded area indicates boundary cells 
to which boundary conditions need to be supplied. The unshaded area, where both i > i m [ n and j > j m i n , is 
a portion of physical domain where a solution is sought. The dots represent the positions of cell-centered 
variables at half-integer locations. 
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'max 



'min 



"mm 



"max 



Fig. 31. — The four corners of the staggered mesh. As in Figure 30, shaded areas indicate boundary 
zones for which boundary conditions need to be supplied. Unshaded areas represent portions of the physical 
domain, where both / m ; n < i < / max +i and j m j n < j < j max +i- The dots represent the locations of cell-centered 
variables at half-integer locations. 
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For any vector quantity a, the procedure is analogous. In the x\ direction, for the x\ component of vector a, 
we use at left-hand boundary, 

[ ai L n -lJ+(l/2) = [°ilw+(l/2) = [ a l] Imm +lJ+(l/2) . ( K3 ) 

where we note that [ci] (min -2j is undefined. At the right-hand boundary we have 

[ ai kax+2J+(l/2) = [ a lL' rallx +l,;+(l/2) = [ a lL raax J+(l/2) • ( K4 ) 

We also need to specify how 02 behaves at the both right- and left-hand boundaries. At the left-hand 
boundary, 

[Vdinto-Vmj = [ a 2]i rain -(l/2) J = [ <J 2]/ min +(l/2) > ; , (K5) 

and, at the right, 

[ ff 2]w+(5/2) J = [ a 2] !mllx +(3/2),; = [°2] imax+ (l/2) J ■ ( K 6) 

2. Periodic Boundary Conditions. Next, we apply periodic boundary conditions in the x\ direction. For 
a scalar, we have at the left-hand boundary, 

[H mln -(3/2),;+(l/2) = M, max -(l/2)J+(l/2) ( K7 ) 

M !min -(l/2)j+(l/2) = M ;max +(l/2)J+(l/2) > ( K8 ) 

while at the right, 

M, raax +(3/2)J+(l/2) = M, min +(l/2)J+(l/2) (K9) 

Mi max +(5/2),y+(l/2) = M« rain +(3/2)J+(l/2) • ( K1 °) 

For vectors, the procedure is again analogous. In the x\ direction, for <5\ at the left-hand boundary, 

[ ai L„-U+(l/2) = [ a l]i max -lj'+(l/2) » ( K11 ) 

[a^j+ii/i) = [ a iL max j+(i/2) ' ( K12 ) 



and at the right, 



For G2 at the left-hand boundary, 



and at the right, 



[ (7 l]w+lj+(l/2) - [ <J l]i mta +l,;+(l/2) > ( K13 ) 

[ (7 l] ; - m ax+2,y+(l/2) = [ a l]i mm +2J+(l/2) • ( K14 ) 

[^^-(3/2) j = [ a 2] !raax -(l/ 2 )J . (K15) 

[0-2] /mta _(l/2)J = [°2] ;max+ (l/ 2 )J , (K16) 

[ a 2]w+(3/2),; = [ <J 2] imm+ (i /2 ), i , (K17) 

[ <7 2] imax +(5/2)J = [ <J 2]i mta+ (3/2)J • ( K18 ) 
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3. Reflecting Boundary Conditions. For scalars, applying reflecting boundary conditions in the x\ 
direction, at the left boundary gives 

Mu-(l/2)j = Wy.+(l/2)J. (K19) 

and at the right-hand boundary, 

[^]i max +(3/2)J = [VL max +(l/2)J . ( K21 ) 

+ (5/2) J = [V], max -(1/2)J • ( K22 ) 

For vector quantities, reflection symmetry requires that the normal component of a vector vanish at planes 
of symmetry, Thus, for reflecting boundary conditions in the x\ direction, we have the following at the 
left-hand boundary for normal component, Ci : 

( K2 3) 

[ a lL„-lJ+(l/2) = - [ <J l]/ min +lJ+(l/2) • ( K24 ) 

Similarly, at the right-hand boundary, we have 

Hw + u + (i/2)=0, (K25) 

[ a lL raax +2J+(l/2) = " [ ff lLj+(l/2) • (K26) 

Since the tangential vector component, 02, is not defined at the x\ planes of symmetry, we do not have to 

apply a zero-value requirement at boundaries. Hence, boundary values are reflected in the same way as 
scalars. At the left-hand boundary, 

,y = [°2] Imin+( i/2)J . ( K27 ) 

[ a 2] ;mm _ (3 /2) J = [°2] Imi n+(3/2)J . ( K2§ ) 

while at the right, 

[ <7 2]; max +(3/2)J = [ a 2] !max+( i /2 )j , (K29) 

t a 2] W+(5/2) J = [ a 2] !max -(l/ 2 )J • (K30) 

4. Flux-Conserving Boundary Conditions. This option is especially useful for problems with spherical 
geometry, when applied to the radial coordinate. Rather than keeping a physical quantity (such as density) 
constant across the boundary, this option allows flux of a physical quantity (such as mass flux) to be kept 
constant. At the outer boundary, a constant flux of iff is maintained when 

Mw+(3/2)J = [H max +(1/2)J ( ! J ,max+(V2 | J j , (K31) 
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Mw+(5/2)J - [VL max +(l/2)j ( L! W+ | V2)J ) ( K32 ) 

\^W+(5/2)J/ 

If the same condition is desired at an inner boundary, for a problem where the computational grid does not 
extend all the way to r = 0, we have 

= M Imm+ (i/2), f [ r 1' m '^ V2)J ) , (K33) 

[V] W -(3/2)J = Mi^-WJ ( m^^ ) ■ (K34) 

Finally, it should be noted that the spatial boundary conditions we have just outlined are also applied 
to each spectral component of radiation quantities. Additionally, with spectral quantities, there is an another 
dimension for which boundary conditions must be applied — the energy dimension. In V2D, we always apply 
energy-space boundary conditions that prevent both inflow below the lowest group (the lower boundary is 
typically at zero anyway) and outflow beyond the highest group. In the case of the latter, we always carry a 
sufficient number of groups such that radiation occupancy in the higher groups is always small. Thus, it is 
largely irrelevant how we treat this upper boundary. 



L. Enforcing the Pauli Exclusion Principal for Neutrinos 

As we discuss in §2.3, it is necessary to enforce the Pauli exclusion principle for neutrinos. Neither of 
the operator split portions of equation (6), the radiation diffusion equation (HI) and the neutrino advection 
equation (48), guarantee that the neutrino occupancy of a given energy state will remain less than unity. This 
reflects the fact that equation (6) is semiclassical in the sense that only the collision integral is treated in a 
quantum-mechanical fashion. Therefore, we must enforce the Pauli exclusion principle as a separate step in 
our algorithm. Other authors (Bruenn 1985; Mezzacappa & Bruenn 1993b) have adopted similar strategies. 

The Pauli exclusion principle requires that the constraint (30) be satisfied. In general, when a numer- 
ical solution of equations (HI) or (48) are obtained, the newly calculated value of the neutrino radiation 
energy density E e may not satisfy this constraint. As neutrinos are produced during the collapse of a stellar 
core, the lower energy states in the center of the core become fully populated. The continued hydrodynamic 
compression of neutrinos within a spatial zone, as described by equation (HI) or (48), can result in occu- 
pation numbers greater than unity in a given group. In our finite difference notation, the Pauli constraint 
corresponds to 

«[£ek+(l/2),-+(l/2)J+(l/2) 



(j £ W(l/2)) 



< 1, (LI) 



where a = (hc) 3 /4ng = 9.4523 MeV 4 cm 3 erg 1 for both photons and neutrinos (for which g = 1). 
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In order to counter this unphysical effect, we examine the radiation energy densities each time equa- 
tion (HI) or (48) is solved. If the values of the radiation energy density in a specific group exceeds e 3 /a, 
the excess neutrinos are removed from that group and placed in the next highest energy group or groups 
where phase space is available. After this process is completed for all groups, the matter internal energy is 
corrected to account for any change in the total (integrated over the spectrum) neutrino energy density. A 
new temperature and pressure are then computed for the zone. 

This enforcement algorithm is applied as the final operation in all substeps corresponding to the red 
boxes in Figure 3. 
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